added source code
[unres.git] / source / unres / src_MD / src / md-diff / np / 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       include 'COMMON.CONTROL'
10       include 'COMMON.VAR'
11       include 'COMMON.MD'
12       include 'COMMON.LANGEVIN'
13       include 'COMMON.CHAIN'
14       include 'COMMON.DERIV'
15       include 'COMMON.GEO'
16       include 'COMMON.LOCAL'
17       include 'COMMON.INTERACT'
18       include 'COMMON.IOUNITS'
19       include 'COMMON.NAMES'
20       include 'COMMON.TIME1'
21       include 'COMMON.MUCA'
22       double precision energia(0:n_ene),vcm(3),incr(3)
23       double precision cm(3),L(3)
24       integer ilen,count,rstcount
25       external ilen
26       character*50 tytul
27       integer maxcount_scale /20/
28       common /gucio/ cm, energia
29       double precision stochforcvec(MAXRES6)
30       common /stochcalc/ stochforcvec
31       integer itime
32       logical scale
33       double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
34       double precision vtnp_(maxres6),vtnp_a(maxres6)
35 c
36       scale=.true.
37       icount_scale=0
38       if (lang.eq.1) then
39         call sddir_precalc
40       else if (lang.eq.2 .or. lang.eq.3) then
41         call stochastic_force(stochforcvec)
42       endif
43       itime_scal=0
44       do while (scale)
45         icount_scale=icount_scale+1
46         if (icount_scale.gt.maxcount_scale) then
47           write (iout,*) 
48      &      "ERROR: too many attempts at scaling down the time step. ",
49      &      "amax=",amax,"epdrift=",epdrift,
50      &      "damax=",damax,"edriftmax=",edriftmax,
51      &      "d_time=",d_time
52             stop
53         endif
54 c First step of the velocity Verlet algorithm
55         if (lang.eq.2) then
56           call sd_verlet1
57         else if (lang.eq.3) then
58           call sd_verlet1_ciccotti
59         else if (lang.eq.1) then
60           call sddir_verlet1
61         else if (tnp1) then
62           call tnp1_step1
63         else if (tnp) then
64           call tnp_step1
65         else    
66           if (tnh) then
67
68             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
69  
70             do i=0,2*nres
71              do j=1,3
72               d_t_old(j,i)=d_t_old(j,i)*scale_nh
73              enddo
74             enddo 
75           endif
76           call verlet1
77         endif     
78 c Build the chain from the newly calculated coordinates 
79         call chainbuild_cart
80         if (rattle) call rattle1
81         if (large.and. mod(itime,ntwe).eq.0) then
82           write (iout,*) "Cartesian and internal coordinates: step 1"
83           call cartprint
84           call intout
85           write (iout,*) "Accelerations"
86           do i=0,nres
87             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
88      &      (d_a(j,i+nres),j=1,3)
89           enddo
90           write (iout,*) "Velocities, step 1"
91           do i=0,nres
92             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
93      &      (d_t(j,i+nres),j=1,3)
94           enddo
95         endif
96         tt0 = tcpu()
97 c Calculate energy and forces
98         call zerograd
99         E_old=potE
100         call etotal(energia)
101         potE=energia(0)-energia(20)
102         call cartgrad
103 c Get the new accelerations
104         call lagrangian
105         t_enegrad=t_enegrad+tcpu()-tt0
106 c Determine maximum acceleration and scale down the timestep if needed
107         call max_accel
108         call predict_edrift(epdrift)
109         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
110 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
111           scale=.true.
112           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
113      &      /dlog(2.0d0)+1
114           itime_scal=itime_scal+ifac_time
115 c          fac_time=dmin1(damax/amax,0.5d0)
116           fac_time=0.5d0**ifac_time
117           d_time=d_time*fac_time
118           if (lang.eq.2 .or. lang.eq.3) then 
119 c            write (iout,*) "Calling sd_verlet_setup: 1"
120 c Rescale the stochastic forces and recalculate or restore 
121 c the matrices of tinker integrator
122             if (itime_scal.gt.maxflag_stoch) then
123               if (large) write (iout,'(a,i5,a)') 
124      &         "Calculate matrices for stochastic step;",
125      &         " itime_scal ",itime_scal
126               if (lang.eq.2) then
127                 call sd_verlet_p_setup
128               else
129                 call sd_verlet_ciccotti_setup
130               endif
131               write (iout,'(2a,i3,a,i3,1h.)') 
132      &         "Warning: cannot store matrices for stochastic",
133      &         " integration because the index",itime_scal,
134      &         " is greater than",maxflag_stoch
135               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
136      &         " integration Langevin algorithm for better efficiency."
137             else if (flag_stoch(itime_scal)) then
138               if (large) write (iout,'(a,i5,a,l1)') 
139      &         "Restore matrices for stochastic step; itime_scal ",
140      &         itime_scal," flag ",flag_stoch(itime_scal)
141               do i=1,dimen
142                 do j=1,dimen
143                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
144                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
145                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
146                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
147                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
148                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
149                 enddo
150               enddo
151             else
152               if (large) write (iout,'(2a,i5,a,l1)') 
153      &         "Calculate & store matrices for stochastic step;",
154      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
155               if (lang.eq.2) then
156                 call sd_verlet_p_setup  
157               else
158                 call sd_verlet_ciccotti_setup
159               endif
160               flag_stoch(ifac_time)=.true.
161               do i=1,dimen
162                 do j=1,dimen
163                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
164                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
165                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
166                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
167                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
168                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
169                 enddo
170               enddo
171             endif
172             fac_time=1.0d0/dsqrt(fac_time)
173             do i=1,dimen
174               stochforcvec(i)=fac_time*stochforcvec(i)
175             enddo
176           else if (lang.eq.1) then
177 c Rescale the accelerations due to stochastic forces
178             fac_time=1.0d0/dsqrt(fac_time)
179             do i=1,dimen
180               d_as_work(i)=d_as_work(i)*fac_time
181             enddo
182           endif
183           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')  
184      &      "itime",itime," Timestep scaled down to ",
185      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
186         else 
187 c Second step of the velocity Verlet algorithm
188           if (lang.eq.2) then   
189             call sd_verlet2
190           else if (lang.eq.3) then
191             call sd_verlet2_ciccotti
192           else if (lang.eq.1) then
193             call sddir_verlet2
194           else if (tnp1) then
195             call tnp1_step2
196           else if (tnp) then
197             call tnp_step2
198           else
199             call verlet2
200             if (tnh) then
201               call kinetic(EK)
202               call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
203               do i=0,2*nres
204                do j=1,3
205                 d_t(j,i)=d_t(j,i)*scale_nh
206                enddo
207               enddo 
208             endif
209           endif                     
210           if (rattle) call rattle2
211           totT=totT+d_time
212           if (d_time.ne.d_time0) then
213             d_time=d_time0
214             if (lang.eq.2 .or. lang.eq.3) then
215               if (large) write (iout,'(a)') 
216      &         "Restore original matrices for stochastic step"
217 c              write (iout,*) "Calling sd_verlet_setup: 2"
218 c Restore the matrices of tinker integrator if the time step has been restored
219               do i=1,dimen
220                 do j=1,dimen
221                   pfric_mat(i,j)=pfric0_mat(i,j,0)
222                   afric_mat(i,j)=afric0_mat(i,j,0)
223                   vfric_mat(i,j)=vfric0_mat(i,j,0)
224                   prand_mat(i,j)=prand0_mat(i,j,0)
225                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
226                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
227                 enddo
228               enddo
229             endif
230           endif
231           scale=.false.
232         endif
233       enddo
234 c Calculate the kinetic and the total energy and the kinetic temperature
235       if (tnp .or. tnp1) then 
236        do i=0,2*nres
237         do j=1,3
238           d_t_old(j,i)=d_t(j,i)
239           d_t(j,i)=d_t(j,i)/s_np
240         enddo
241        enddo 
242       endif
243       call kinetic(EK)
244       totE=EK+potE
245 c diagnostics
246 c      call kinetic1(EK1)
247 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
248 c end diagnostics
249 c Couple the system to Berendsen bath if needed
250       if (tbf .and. lang.eq.0) then
251         call verlet_bath
252       endif
253       kinetic_T=2.0d0/(dimen*Rb)*EK
254 c Backup the coordinates, velocities, and accelerations
255       do i=0,2*nres
256         do j=1,3
257           dc_old(j,i)=dc(j,i)
258           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
259           d_a_old(j,i)=d_a(j,i)
260         enddo
261       enddo 
262 c      if (mod(itime,ntwe).eq.0 .and. large) then
263       if (mod(itime,ntwe).eq.0) then
264
265        if(tnp .or. tnp1) then
266         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
267         H=(HNose1-H0)*s_np
268 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
269 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
270         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
271        endif
272
273        if(tnh) then
274         HNose1=Hnose_nh(EK,potE)
275         H=HNose1-H0
276         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
277        endif
278
279        if (large) then
280         itnp=0
281         do j=1,3
282          itnp=itnp+1
283          vtnp(itnp)=d_t(j,0)
284         enddo
285         do i=nnt,nct-1  
286          do j=1,3    
287           itnp=itnp+1
288           vtnp(itnp)=d_t(j,i)
289          enddo
290         enddo
291         do i=nnt,nct
292          if (itype(i).ne.10) then
293           inres=i+nres
294           do j=1,3  
295            itnp=itnp+1  
296            vtnp(itnp)=d_t(j,inres)
297           enddo
298          endif      
299         enddo 
300
301 c Transform velocities from UNRES coordinate space to cartesian and Gvec
302 c eigenvector space
303
304         do i=1,dimen
305           vtnp_(i)=0.0d0
306           vtnp_a(i)=0.0d0
307           do j=1,dimen
308             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
309             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
310           enddo
311           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
312         enddo
313
314         do i=1,dimen
315          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
316         enddo
317
318        endif
319       endif
320       return
321       end