Statistics
| Revision:

root / src / egrad_mopac.f90 @ 10

History | View | Annotate | Download (4.9 kB)

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
       FPBC, Lat_a, Lat_b, Lat_c, IPer
8
  use Io_module
9

    
10
  IMPLICIT NONE
11

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

    
19
  ! ======================================================================
20

    
21
  character(LCHARS) :: LINE
22

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

    
26
  REAL(KREAL) :: Pi
27

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

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

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

    
36
  ! ======================================================================
37

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

    
45
  ! ======================================================================
46

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

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

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

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

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

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

    
70
  ! we create the input file
71

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

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

    
87
  DO I=1,Nat
88
     If (renum) THEN
89
        Iat=Order(I)
90
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:)
91
     ELSE
92
        Iat=OrderInv(I)
93
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:)
94
     END IF
95
  END DO
96
! If we have periodic boundaries, we have to print them
97
  IF (FPBC) THEN
98
     IF (IPER>=1) THEN
99
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_a(1:3)
100
     END IF
101
     IF (IPER>=2) THEN
102
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_b(1:3)
103
     END IF
104
     IF (IPER>=3) THEN
105
        WRITE(IOTMP,'(1X," Tv ",3(1X,F15.8))') Lat_c(1:3)
106
     END IF
107
  END IF
108

    
109
     WRITE(IOTMP,'(1X,A)') Trim(Mopac_EndGeom)
110

    
111
  Current => Mopac_End
112
  DO WHILE (ASSOCIATED(Current%next))
113
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
114
     Current => current%next
115
  END DO
116

    
117
  WRITE(IOTMP,*) 
118
  CLOSE(IOTMP)
119

    
120
  IF (debug) THEN
121
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
122
     call system('cat ' // Trim(FileIn))
123
  END IF
124

    
125

    
126
  call system(RunCommand)
127

    
128
  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'
129

    
130
  ! Whatever the coordinate system, we read the cartesian gradient
131
  ! and we convert it to Zmat or other if needed
132

    
133
  OPEN(IOTMP,FILE=FileOut, STATUS='old')
134
  ! We first search for the forces
135
  READ(IOTMP,'(A)') LINE
136

    
137
  DO WHILE (INDEX(LINE,'FINAL  POINT  AND  DERIVATIVES')==0)
138
     IF (INDEX(Line,'FINAL HEAT OF FORMATION').NE.0)  THEN
139
        Itmp1=Index(Line,"=")+1
140
        ITmp2=Index(Line,"KCAL")-1
141
        Read(Line(Itmp1:Itmp2),*) E
142
        E=E/au2kcal
143
     END IF
144
     READ(IOTMP,'(A)',END=999) LINE
145
  END DO
146

    
147
  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
148
  READ(IOTMP,'(A)') LINE
149
  READ(IOTMP,'(A)') LINE
150
  DO I=1,Nat
151
     Iat=I
152
     IF (renum) Iat=Order(I)
153
     DO J=1,3
154
        READ(IOTMP,'(A)') Line
155
! Gradient is on the 7th column
156
        DO K=1,6
157
           Line=AdjustL(Line)
158
           Idx=Index(Line," ")
159
           Line=AdjustL(Line(Idx+1:))
160
        END DO
161
        READ(Line,*) GradCart(3*Iat-3+j)
162
     END DO
163
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
164
  END DO
165

    
166
! Mopac gives the Forces in kcal/angstrom, so we convert it into the 
167
! gradient in ua/Angstrom 
168
  gradCart=gradCart/au2kcal
169

    
170
!  if (debug.AND.(COORD.EQ.'ZMAT')) THEN
171
!     DO WHILE (INDEX(LINE,'Internal Coordinate Forces')==0)
172
!        READ(IOTMP,'(A)') LINE
173
!!        WRITE(*,*) 'Pas bon:' // TRIM(LINE)
174
!     END DO
175
!     WRITE(*,*) TRIM(LINE)
176
!     DO I=1,Nat+2
177
!        READ(IOTMP,'(A)') LINE
178
!        WRITE(*,*) TRIM(LINE)
179
!     END DO
180
!  END IF
181

    
182
  CLOSE(IOTMP)
183

    
184
  if (debug) WRITE(*,*) '================  Egrad_Mopac Over ===================='
185

    
186
  RETURN
187

    
188
999 CONTINUE
189
  WRITE(*,*) 'We should not be here !!!!'
190
  STOP
191
  ! ======================================================================
192
end subroutine egrad_mopac