Statistiques
| Révision :

root / src / egrad_mopac.f90 @ 8

Historique | Voir | Annoter | Télécharger (4,54 ko)

1
subroutine egrad_mopac(e,geomcart,gradcart)
2

    
3
  ! This routines calculates the energy and the gradient of 
4
  ! a molecule, using mopac 
5

    
6
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, ProgExe
7
  use Io_module
8

    
9
  IMPLICIT NONE
10

    
11
  ! Energy (calculated if F300K=.F., else estimated)
12
  REAL(KREAL), INTENT (OUT) :: e
13
  ! Nb degree of freedom
14
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
15
  ! Gradient calculated at Geom geometry
16
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)
17

    
18
  ! ======================================================================
19

    
20
  character(LCHARS) :: LINE
21

    
22
  logical           :: debug
23
  LOGICAL, SAVE :: first=.true.
24

    
25
  REAL(KREAL) :: Pi
26

    
27
  INTEGER(KINT) :: iat, i, j, k, n3at
28
  INTEGER(KINT) :: ITmp1,ITmp2, Idx
29

    
30
  CHARACTER(132) :: FileIn, FileOut, FileChk
31

    
32
  CHARACTER(3) :: SepIn=' < ', SepOut=' > '
33
  CHARACTER(VLCHARS), SAVE ::  RunCommand
34

    
35
  ! ======================================================================
36

    
37
  INTERFACE
38
     function valid(string) result (isValid)
39
       CHARACTER(*), intent(in) :: string
40
       logical                  :: isValid
41
     END function VALID
42
  END INTERFACE
43

    
44
  ! ======================================================================
45

    
46
  Pi=dacos(-1.0d0)
47
  n3at=3*nat
48

    
49
  debug=valid('EGRAD')
50
  if (debug) WRITE(*,*) '================ Entering Egrad_Mopac ===================='
51

    
52
  if (first) THEN
53
     First=.FALSE.
54
     SepIn='   '
55
     SepOut='   '
56
  END IF
57
 
58
  gradcart=0.
59
  E=0.
60

    
61
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)
62
  FileOut=Trim(CalcName) // "." // Trim(OSuffix)
63
  FileChk=Trim(CalcName) // ".arc"
64

    
65
  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) 
66

    
67
  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)
68

    
69
  ! we create the input file
70

    
71
  OPEN(IOTMP,File=FileIn)
72
  WRITE(IOTMP,'(A)') '* All comments are put at the begining of the file'
73
  Current => Mopac_comment
74
  DO WHILE (ASSOCIATED(Current%next))
75
     WRITE(IOTMP,'(A)') Trim(current%line)
76
     Current => current%next
77
  END DO
78

    
79
  Current => Mopac_root
80
  DO WHILE (ASSOCIATED(Current%next))
81
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
82
     WRITE(*,'(1X,A)') Trim(current%line)
83
     Current => current%next
84
  END DO
85

    
86
  DO I=1,Nat
87
     If (renum) THEN
88
        Iat=Order(I)
89
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:)
90
     ELSE
91
        Iat=OrderInv(I)
92
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:)
93
     END IF
94
  END DO
95

    
96
     WRITE(IOTMP,'(1X,A)') Trim(Mopac_EndGeom)
97

    
98
  Current => Mopac_End
99
  DO WHILE (ASSOCIATED(Current%next))
100
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
101
     Current => current%next
102
  END DO
103

    
104
  WRITE(IOTMP,*) 
105
  CLOSE(IOTMP)
106

    
107
  IF (debug) THEN
108
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
109
     call system('cat ' // Trim(FileIn))
110
  END IF
111

    
112

    
113
  call system(RunCommand)
114

    
115
  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'
116

    
117
  ! Whatever the coordinate system, we read the cartesian gradient
118
  ! and we convert it to Zmat or other if needed
119

    
120
  OPEN(IOTMP,FILE=FileOut, STATUS='old')
121
  ! We first search for the forces
122
  READ(IOTMP,'(A)') LINE
123

    
124
  DO WHILE (INDEX(LINE,'FINAL  POINT  AND  DERIVATIVES')==0)
125
     IF (INDEX(Line,'FINAL HEAT OF FORMATION').NE.0)  THEN
126
        Itmp1=Index(Line,"=")+1
127
        ITmp2=Index(Line,"KCAL")-1
128
        Read(Line(Itmp1:Itmp2),*) E
129
        E=E/au2kcal
130
     END IF
131
     READ(IOTMP,'(A)',END=999) LINE
132
  END DO
133

    
134
  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
135
  READ(IOTMP,'(A)') LINE
136
  READ(IOTMP,'(A)') LINE
137
  DO I=1,Nat
138
     Iat=I
139
     IF (renum) Iat=Order(I)
140
     DO J=1,3
141
        READ(IOTMP,'(A)') Line
142
! Gradient is on the 7th column
143
        DO K=1,6
144
           Line=AdjustL(Line)
145
           Idx=Index(Line," ")
146
           Line=AdjustL(Line(Idx+1:))
147
        END DO
148
        READ(Line,*) GradCart(3*Iat-3+j)
149
     END DO
150
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
151
  END DO
152

    
153
! Mopac gives the Forces in kcal/angstrom, so we convert it into the 
154
! gradient in ua/Angstrom 
155
  gradCart=gradCart/au2kcal
156

    
157
!  if (debug.AND.(COORD.EQ.'ZMAT')) THEN
158
!     DO WHILE (INDEX(LINE,'Internal Coordinate Forces')==0)
159
!        READ(IOTMP,'(A)') LINE
160
!!        WRITE(*,*) 'Pas bon:' // TRIM(LINE)
161
!     END DO
162
!     WRITE(*,*) TRIM(LINE)
163
!     DO I=1,Nat+2
164
!        READ(IOTMP,'(A)') LINE
165
!        WRITE(*,*) TRIM(LINE)
166
!     END DO
167
!  END IF
168

    
169
  CLOSE(IOTMP)
170

    
171
  if (debug) WRITE(*,*) '================  Egrad_Mopac Over ===================='
172

    
173
  RETURN
174

    
175
999 CONTINUE
176
  WRITE(*,*) 'We should not be here !!!!'
177
  STOP
178
  ! ======================================================================
179
end subroutine egrad_mopac