1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
|
**********************************************************************
* *
* REAXFF Reactive force field program *
* *
* Developed and written by Adri van Duin, duin@wag.caltech.edu *
* *
* Copyright (c) 2001-2010 California Institute of Technology *
* *
* This is an open-source program. Feel free to modify its *
* contents. Please keep me informed of any useful modification *
* or addition that you made. Please do not distribute this *
* program to others; if people are interested in obtaining *
* a copy of this program let them contact me first. *
* *
**********************************************************************
**********************************************************************
* *
* REAXFF Reactive force field program *
* *
* Developed and written by Adri van Duin, duin@wag.caltech.edu *
* *
* Copyright (c) 2001-2010 California Institute of Technology *
* *
* This is an open-source program. Feel free to modify its *
* contents. Please keep me informed of any useful modification *
* or addition that you made. Please do not distribute this *
* program to others; if people are interested in obtaining *
* a copy of this program let them contact me first. *
* *
**********************************************************************
**********************************************************************
subroutine encalc
**********************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkcha.blk"
#include "cbkconst.blk"
#include "cbkd.blk"
#include "cbkdcell.blk"
#include "cbkenergies.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "small.blk"
**********************************************************************
* *
* Calculate energy and first derivatives *
* *
**********************************************************************
c$$$ if (ndebug.eq.1) then
c$$$C open (65,file='fort.65',status='unknown',access='append')
c$$$ write (65,*) 'In encalc'
c$$$ call timer(65)
c$$$ close (65)
c$$$ end if
estrc=0.0
do i1=1,na
do i2=1,3
d(i2,i1)=0.0
estrain(i1)=0.0
end do
end do
eb=zero
ea=zero
elp=zero
emol=zero
ev=zero
ehb=zero
ecoa=zero
epen=zero
et=zero
eco=zero
eres=zero
eradbo=zero
efi=zero
if(Lvirial.eq.1) then
do k1 = 1,6
virial(k1) = zero
end do
endif
if (Latomvirial.eq.1) then
do i1=1,na
do i2=1,6
atomvirial(i2,i1)=0.0
end do
end do
endif
call boncor
call lonpar
call covbon
call ovcor
call srtang !Determine valency angles
call srttor !Determine torsion angles
* call srtoop !Determine out of plane angles
call srthb !Determine hydrogen bonds
call calval
call valang
* call oopang
call torang
call hbond
!print *, 'called hbond'
!print *, nchaud
call nonbon
call efield
call restraint
c
c Use this to print out fort.73-style energies
c It only works correctly in serial mode
c
c write (6,'(i8,1x,14(f21.10,1x))')mdstep+nprevrun,eb,ea,elp,
c $emol,ev+epen,ecoa,ehb,et,eco,ew,ep,ech,efi
estrc=eb+ea+elp+ev+ecoa+emol+epen+et+ehb+eco+ew+ep+ncha2*ech+efi
if (estrc.gt.zero) return
if (estrc.le.zero) then
goto 10
else
write (*,*)mdstep
write (92,*)eb,ea,elp,ev,ecoa,emol,epen,eoop,et,eco,ew,
$ep,ech,eres,eradbo
stop 'Energy not a number'
end if
10 continue
return
end
**********************************************************************
**********************************************************************
subroutine reaxinit
**********************************************************************
#include "cbka.blk"
#include "cbkatomcoord.blk"
#include "cbkc.blk"
#include "cbkcha.blk"
#include "cbkconst.blk"
#include "cbkdcell.blk"
#include "cbkdistan.blk"
#include "cbkenergies.blk"
#include "cbkia.blk"
#include "cbkimove.blk"
#include "cbkinit.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
**********************************************************************
c$$$ if (ndebug.eq.1) then
c$$$C open (65,file='fort.65',status='unknown',access='append')
c$$$ write (65,*) 'In init'
c$$$ call timer(65)
c$$$ close (65)
c$$$ end if
**********************************************************************
* *
* Initialize variables *
* *
**********************************************************************
convmd=4.184*1.0d26
pi=3.14159265
avognr=6.0221367d23
rdndgr=180.0/pi
dgrrdn=1.0/rdndgr
rgasc=8.314510
caljou=4.184
xjouca=1.0/caljou
ech=zero
zero=0.0
one=1.0
two=2.0
three=3.0
half=one/two
nzero=0
none=1
ntwo=2
nthree=3
invt=0
ndata2=0
iheatf=0
nradcount=0
itemp=1
xinh=zero
ifieldx=0
ifieldy=0
ifieldz=0
mdstep=0
kx=0
ky=0
kz=0
nit=0
nbon=0
angle(1)=90.0
angle(2)=90.0
angle(3)=90.0
axiss(1)=zero
axiss(2)=zero
axiss(3)=zero
do i1=1,nat
id(i1,1)=0
id(i1,2)=0
id(i1,3)=0
end do
icgeo=0
sumhe=zero
ustime=zero
systime=zero
vpmax=zero
vpmin=zero
dseed=0
iagain=0
do i1=1,nat
do i2=1,mbond+3
ia(i1,i2)=0
iag(i1,i2)=0
end do
end do
ioldchg=0
na=0
nrestra=0
nrestrav=0
nrestrat=0
nrestram=0
tset=tsetor
tm11=axis(1)
tm21=zero
tm31=zero
tm22=axis(2)
tm32=zero
tm33=axis(3)
qruid='NORMAL RUN'
c$$$ do i1=1,nat
c$$$ imove(i1)=1
c$$$ end do
**********************************************************************
* *
* Write file headers *
* *
**********************************************************************
Cc open (71,file='fort.71',status='unknown',access='append')
c write (71,10)
c close (71)
Cc open (73,file='fort.73',status='unknown',access='append')
c write (73,20)
c close (73)
c if (ntrc.gt.0) then
Cc open (75,file='fort.75',status='unknown',access='append')
c write (75,30)
c close (75)
c end if
c if (nmethod.eq.4) then
Cc open (59,file='fort.59',status='unknown',access='append')
c write (59,40)
c close (59)
c end if
return
**********************************************************************
* *
* Format part *
* *
**********************************************************************
10 format (' Iter. Nmol Epot Ekin Etot ',
$' T(K) Eaver(block) Eaver(total) Taver Tmax ',
$' Pres(MPa) sdev(Epot) sdev(Eaver) Tset Timestep',
$' RMSG Totaltime')
20 format (' Iter. Ebond Eatom Elp Emol',
$' Eval Ecoa Ehbo Etors Econj',
$' Evdw Ecoul Echarge Efield')
30 format (' Iter. Tsys Tzone1 Tset1 Tzone2 Tset2')
40 format (' Iter. a b c px',
$'(MPa) py(MPa) pz(MPa) pset(MPa) Volume ')
end
**********************************************************************
************************************************************************
c subroutine timer(nunit)
************************************************************************
c#include "cbka.blk"
c#include "cbkinit.blk"
c real timear
c real tarray(2)
c#ifdef _IBM
c call dtime_(tarray,timear)
c#else
c call dtime(tarray,timear)
c#endif
c ustime=ustime+tarray(1)
c systime=systime+tarray(2)
c write (nunit,100)ustime,systime,ustime+systime
c return
c 100 format ('User time:',f20.4,' System time:',f20.4,
c $' Total time:',f20.4)
c end
************************************************************************
************************************************************************
|
