added source code
[unres.git] / source / unres / src_MD / md-diff / mts / velverlet_step.f
1 c-------------------------------------------------------------------------------
2       subroutine velverlet_step(itime)
3 c-------------------------------------------------------------------------------
4 c  Perform a single velocity Verlet step; the time step can be rescaled if 
5 c  increments in accelerations exceed the threshold
6 c-------------------------------------------------------------------------------
7       implicit real*8 (a-h,o-z)
8       include 'DIMENSIONS'
9 #ifdef MPI
10       include 'mpif.h'
11       integer ierror,ierrcode
12 #endif
13       include 'COMMON.SETUP'
14       include 'COMMON.CONTROL'
15       include 'COMMON.VAR'
16       include 'COMMON.MD'
17 #ifndef LANG0
18       include 'COMMON.LANGEVIN'
19 #else
20       include 'COMMON.LANGEVIN.lang0'
21 #endif
22       include 'COMMON.CHAIN'
23       include 'COMMON.DERIV'
24       include 'COMMON.GEO'
25       include 'COMMON.LOCAL'
26       include 'COMMON.INTERACT'
27       include 'COMMON.IOUNITS'
28       include 'COMMON.NAMES'
29       include 'COMMON.TIME1'
30       include 'COMMON.MUCA'
31       double precision vcm(3),incr(3)
32       double precision cm(3),L(3)
33       integer ilen,count,rstcount
34       external ilen
35       character*50 tytul
36       integer maxcount_scale /20/
37       common /gucio/ cm
38       double precision stochforcvec(MAXRES6)
39       common /stochcalc/ stochforcvec
40       integer itime
41       logical scale
42 c
43       scale=.true.
44       icount_scale=0
45       if (lang.eq.1) then
46         call sddir_precalc
47       else if (lang.eq.2 .or. lang.eq.3) then
48 #ifndef LANG0
49         call stochastic_force(stochforcvec)
50 #else
51         write (iout,*) 
52      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
53 #ifdef MPI
54         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
55 #endif
56         stop
57 #endif
58       endif
59       itime_scal=0
60       do while (scale)
61         icount_scale=icount_scale+1
62         if (icount_scale.gt.maxcount_scale) then
63           write (iout,*) 
64      &      "ERROR: too many attempts at scaling down the time step. ",
65      &      "amax=",amax,"epdrift=",epdrift,
66      &      "damax=",damax,"edriftmax=",edriftmax,
67      &      "d_time=",d_time
68           call flush(iout)
69 #ifdef MPI
70           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
71 #endif
72           stop
73         endif
74 c First step of the velocity Verlet algorithm
75         if (lang.eq.2) then
76 #ifndef LANG0
77           call sd_verlet1
78 #endif
79         else if (lang.eq.3) then
80 #ifndef LANG0
81           call sd_verlet1_ciccotti
82 #endif
83         else if (lang.eq.1) then
84           call sddir_verlet1
85         else
86           call verlet1
87         endif     
88 c Build the chain from the newly calculated coordinates 
89         call chainbuild_cart
90         if (rattle) call rattle1
91         if (ntwe.ne.0) then
92         if (large.and. mod(itime,ntwe).eq.0) then
93           write (iout,*) "Cartesian and internal coordinates: step 1"
94           call cartprint
95           call intout
96           write (iout,*) "dC"
97           do i=0,nres
98             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
99      &      (dc(j,i+nres),j=1,3)
100           enddo
101           write (iout,*) "Accelerations"
102           do i=0,nres
103             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
104      &      (d_a(j,i+nres),j=1,3)
105           enddo
106           write (iout,*) "Velocities, step 1"
107           do i=0,nres
108             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
109      &      (d_t(j,i+nres),j=1,3)
110           enddo
111         endif
112         endif
113 #ifdef MPI
114         tt0 = MPI_Wtime()
115 #else
116         tt0 = tcpu()
117 #endif
118 c Calculate energy and forces
119         call zerograd
120         call etotal(potEcomp)
121 #ifdef TIMING_ENE
122 #ifdef MPI
123         t_etotal=t_etotal+MPI_Wtime()-tt0
124 #else
125         t_etotal=t_etotal+tcpu()-tt0
126 #endif
127 #endif
128         potE=potEcomp(0)-potEcomp(20)
129         call cartgrad
130 c Get the new accelerations
131         call lagrangian
132 #ifdef MPI
133         t_enegrad=t_enegrad+MPI_Wtime()-tt0
134 #else
135         t_enegrad=t_enegrad+tcpu()-tt0
136 #endif
137 c Determine maximum acceleration and scale down the timestep if needed
138         call max_accel
139         amax=amax/(itime_scal+1)**2
140         call predict_edrift(epdrift)
141         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
142 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
143           scale=.true.
144           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
145      &      /dlog(2.0d0)+1
146           itime_scal=itime_scal+ifac_time
147 c          fac_time=dmin1(damax/amax,0.5d0)
148           fac_time=0.5d0**ifac_time
149           d_time=d_time*fac_time
150           if (lang.eq.2 .or. lang.eq.3) then 
151 #ifndef LANG0
152 c            write (iout,*) "Calling sd_verlet_setup: 1"
153 c Rescale the stochastic forces and recalculate or restore 
154 c the matrices of tinker integrator
155             if (itime_scal.gt.maxflag_stoch) then
156               if (large) write (iout,'(a,i5,a)') 
157      &         "Calculate matrices for stochastic step;",
158      &         " itime_scal ",itime_scal
159               if (lang.eq.2) then
160                 call sd_verlet_p_setup
161               else
162                 call sd_verlet_ciccotti_setup
163               endif
164               write (iout,'(2a,i3,a,i3,1h.)') 
165      &         "Warning: cannot store matrices for stochastic",
166      &         " integration because the index",itime_scal,
167      &         " is greater than",maxflag_stoch
168               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
169      &         " integration Langevin algorithm for better efficiency."
170             else if (flag_stoch(itime_scal)) then
171               if (large) write (iout,'(a,i5,a,l1)') 
172      &         "Restore matrices for stochastic step; itime_scal ",
173      &         itime_scal," flag ",flag_stoch(itime_scal)
174               do i=1,dimen
175                 do j=1,dimen
176                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
177                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
178                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
179                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
180                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
181                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
182                 enddo
183               enddo
184             else
185               if (large) write (iout,'(2a,i5,a,l1)') 
186      &         "Calculate & store matrices for stochastic step;",
187      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
188               if (lang.eq.2) then
189                 call sd_verlet_p_setup  
190               else
191                 call sd_verlet_ciccotti_setup
192               endif
193               flag_stoch(ifac_time)=.true.
194               do i=1,dimen
195                 do j=1,dimen
196                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
197                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
198                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
199                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
200                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
201                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
202                 enddo
203               enddo
204             endif
205             fac_time=1.0d0/dsqrt(fac_time)
206             do i=1,dimen
207               stochforcvec(i)=fac_time*stochforcvec(i)
208             enddo
209 #endif
210           else if (lang.eq.1) then
211 c Rescale the accelerations due to stochastic forces
212             fac_time=1.0d0/dsqrt(fac_time)
213             do i=1,dimen
214               d_as_work(i)=d_as_work(i)*fac_time
215             enddo
216           endif
217           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')  
218      &      "itime",itime," Timestep scaled down to ",
219      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
220         else 
221 c Second step of the velocity Verlet algorithm
222           if (lang.eq.2) then   
223 #ifndef LANG0
224             call sd_verlet2
225 #endif
226           else if (lang.eq.3) then
227 #ifndef LANG0
228             call sd_verlet2_ciccotti
229 #endif
230           else if (lang.eq.1) then
231             call sddir_verlet2
232           else
233             call verlet2
234           endif                     
235           if (rattle) call rattle2
236           totT=totT+d_time
237           if (d_time.ne.d_time0) then
238             d_time=d_time0
239 #ifndef   LANG0
240             if (lang.eq.2 .or. lang.eq.3) then
241               if (large) write (iout,'(a)') 
242      &         "Restore original matrices for stochastic step"
243 c              write (iout,*) "Calling sd_verlet_setup: 2"
244 c Restore the matrices of tinker integrator if the time step has been restored
245               do i=1,dimen
246                 do j=1,dimen
247                   pfric_mat(i,j)=pfric0_mat(i,j,0)
248                   afric_mat(i,j)=afric0_mat(i,j,0)
249                   vfric_mat(i,j)=vfric0_mat(i,j,0)
250                   prand_mat(i,j)=prand0_mat(i,j,0)
251                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
252                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
253                 enddo
254               enddo
255             endif
256 #endif
257           endif
258           scale=.false.
259         endif
260       enddo
261 c Calculate the kinetic and the total energy and the kinetic temperature
262       call kinetic(EK)
263       totE=EK+potE
264 c diagnostics
265 c      call kinetic1(EK1)
266 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
267 c end diagnostics
268 c Couple the system to Berendsen bath if needed
269       if (tbf .and. lang.eq.0) then
270         call verlet_bath
271       endif
272       kinetic_T=2.0d0/(dimen3*Rb)*EK
273 c Backup the coordinates, velocities, and accelerations
274       do i=0,2*nres
275         do j=1,3
276           dc_old(j,i)=dc(j,i)
277           d_t_old(j,i)=d_t(j,i)
278           d_a_old(j,i)=d_a(j,i)
279         enddo
280       enddo 
281       if (ntwe.ne.0) then
282       if (mod(itime,ntwe).eq.0 .and. large) then
283         write (iout,*) "Velocities, step 2"
284         do i=0,nres
285           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
286      &      (d_t(j,i+nres),j=1,3)
287         enddo
288       endif
289       endif
290       return
291       end