Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / source / unres / src_MD / MD_A-MTS.F
1       subroutine MD
2 c------------------------------------------------
3 c  The driver for molecular dynamics subroutines
4 c------------------------------------------------
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7 #ifdef MPI
8       include "mpif.h"
9       integer IERROR,ERRCODE
10 #endif
11       include 'COMMON.SETUP'
12       include 'COMMON.CONTROL'
13       include 'COMMON.VAR'
14       include 'COMMON.MD'
15 #ifndef LANG0
16       include 'COMMON.LANGEVIN'
17 #else
18       include 'COMMON.LANGEVIN.lang0'
19 #endif
20       include 'COMMON.CHAIN'
21       include 'COMMON.DERIV'
22       include 'COMMON.GEO'
23       include 'COMMON.LOCAL'
24       include 'COMMON.INTERACT'
25       include 'COMMON.IOUNITS'
26       include 'COMMON.NAMES'
27       include 'COMMON.TIME1'
28       double precision cm(3),L(3),vcm(3)
29 #ifdef VOUT
30       double precision v_work(maxres6),v_transf(maxres6)
31 #endif
32       integer ilen,rstcount
33       external ilen
34       character*50 tytul
35       common /gucio/ cm
36       integer itime
37 c
38 #ifdef MPI
39       if (ilen(tmpdir).gt.0)
40      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
41      &        //liczba(:ilen(liczba))//'.rst')
42 #else
43       if (ilen(tmpdir).gt.0)
44      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
45 #endif
46       t_MDsetup=0.0d0
47       t_langsetup=0.0d0
48       t_MD=0.0d0
49       t_enegrad=0.0d0
50       t_sdsetup=0.0d0
51       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
52 #ifdef MPI
53       tt0=MPI_Wtime()
54 #else
55       tt0 = tcpu()
56 #endif
57 c Determine the inverse of the inertia matrix.
58       call setup_MD_matrices
59 c Initialize MD
60       call init_MD
61 #ifdef MPI
62       t_MDsetup = MPI_Wtime()-tt0
63 #else
64       t_MDsetup = tcpu()-tt0
65 #endif
66       rstcount=0 
67 c   Entering the MD loop       
68 #ifdef MPI
69       tt0 = MPI_Wtime()
70 #else
71       tt0 = tcpu()
72 #endif
73       if (lang.eq.2 .or. lang.eq.3) then
74 #ifndef   LANG0
75         call setup_fricmat
76         if (lang.eq.2) then
77           call sd_verlet_p_setup        
78         else
79           call sd_verlet_ciccotti_setup
80         endif
81         do i=1,dimen3
82           do j=1,dimen3
83             pfric0_mat(i,j,0)=pfric_mat(i,j)
84             afric0_mat(i,j,0)=afric_mat(i,j)
85             vfric0_mat(i,j,0)=vfric_mat(i,j)
86             prand0_mat(i,j,0)=prand_mat(i,j)
87             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
88             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
89           enddo
90         enddo
91         flag_stoch(0)=.true.
92         do i=1,maxflag_stoch
93           flag_stoch(i)=.false.
94         enddo  
95 #else
96         write (iout,*) 
97      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
98 #ifdef MPI
99         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
100 #endif
101         stop
102 #endif
103       else if (lang.eq.1 .or. lang.eq.4) then
104         call setup_fricmat
105       endif
106 #ifdef MPI
107       t_langsetup=MPI_Wtime()-tt0
108       tt0=MPI_Wtime()
109 #else
110       t_langsetup=tcpu()-tt0
111       tt0=tcpu()
112 #endif
113       do itime=1,n_timestep
114         rstcount=rstcount+1
115         if (lang.gt.0 .and. surfarea .and. 
116      &      mod(itime,reset_fricmat).eq.0) then
117           if (lang.eq.2 .or. lang.eq.3) then
118 #ifndef LANG0
119             call setup_fricmat
120             if (lang.eq.2) then
121               call sd_verlet_p_setup
122             else
123               call sd_verlet_ciccotti_setup
124             endif
125             do i=1,dimen3
126               do j=1,dimen3
127                 pfric0_mat(i,j,0)=pfric_mat(i,j)
128                 afric0_mat(i,j,0)=afric_mat(i,j)
129                 vfric0_mat(i,j,0)=vfric_mat(i,j)
130                 prand0_mat(i,j,0)=prand_mat(i,j)
131                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
132                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
133               enddo
134             enddo
135             flag_stoch(0)=.true.
136             do i=1,maxflag_stoch
137               flag_stoch(i)=.false.
138             enddo   
139 #endif
140           else if (lang.eq.1 .or. lang.eq.4) then
141             call setup_fricmat
142           endif
143           write (iout,'(a,i10)') 
144      &      "Friction matrix reset based on surface area, itime",itime
145         endif
146         if (reset_vel .and. tbf .and. lang.eq.0 
147      &      .and. mod(itime,count_reset_vel).eq.0) then
148           call random_vel
149           write(iout,'(a,f20.2)') 
150      &     "Velocities reset to random values, time",totT       
151           do i=0,2*nres
152             do j=1,3
153               d_t_old(j,i)=d_t(j,i)
154             enddo
155           enddo
156         endif
157         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
158           call inertia_tensor  
159           call vcm_vel(vcm)
160           do j=1,3
161              d_t(j,0)=d_t(j,0)-vcm(j)
162           enddo
163           call kinetic(EK)
164           kinetic_T=2.0d0/(dimen3*Rb)*EK
165           scalfac=dsqrt(T_bath/kinetic_T)
166           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT       
167           do i=0,2*nres
168             do j=1,3
169               d_t_old(j,i)=scalfac*d_t(j,i)
170             enddo
171           enddo
172         endif  
173         if (lang.ne.4) then
174           if (RESPA) then
175 c Time-reversible RESPA algorithm 
176 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
177             call RESPA_step(itime)
178           else
179 c Variable time step algorithm.
180             call velverlet_step(itime)
181           endif
182         else
183 #ifdef BROWN
184           call brown_step(itime)
185 #else
186           print *,"Brown dynamics not here!"
187 #ifdef MPI
188           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
189 #endif
190           stop
191 #endif
192         endif
193         if (ntwe.ne.0) then
194          if (mod(itime,ntwe).eq.0) call statout(itime)
195 #ifdef VOUT
196         do j=1,3
197           v_work(j)=d_t(j,0)
198         enddo
199         ind=3
200         do i=nnt,nct-1
201           do j=1,3
202             ind=ind+1
203             v_work(ind)=d_t(j,i)
204           enddo
205         enddo
206         do i=nnt,nct
207           if (itype(i).ne.10) then
208             do j=1,3
209               ind=ind+1
210               v_work(ind)=d_t(j,i+nres)
211             enddo
212           endif
213         enddo
214
215         write (66,'(80f10.5)') 
216      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
217         do i=1,ind
218           v_transf(i)=0.0d0
219           do j=1,ind
220             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
221           enddo
222            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
223         enddo
224         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
225 #endif
226         endif
227         if (mod(itime,ntwx).eq.0) then
228           write (tytul,'("time",f8.2)') totT
229           if(mdpdb) then
230              call pdbout(potE,tytul,ipdb)
231           else 
232              call cartout(totT)
233           endif
234         endif
235         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
236            open(irest2,file=rest2name,status='unknown')
237            write(irest2,*) totT,EK,potE,totE,t_bath
238            do i=1,2*nres
239             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
240            enddo
241            do i=1,2*nres
242             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
243            enddo
244           close(irest2)
245           rstcount=0
246         endif 
247       enddo
248 #ifdef MPI
249       t_MD=MPI_Wtime()-tt0
250 #else
251       t_MD=tcpu()-tt0
252 #endif
253       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
254      &  '  Timing  ',
255      & 'MD calculations setup:',t_MDsetup,
256      & 'Energy & gradient evaluation:',t_enegrad,
257      & 'Stochastic MD setup:',t_langsetup,
258      & 'Stochastic MD step setup:',t_sdsetup,
259      & 'MD steps:',t_MD
260       write (iout,'(/28(1h=),a25,27(1h=))') 
261      & '  End of MD calculation  '
262 #ifdef TIMING_ENE
263       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
264      &  " eshort",t_eshort
265       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
266      & " time_fricmatmult",time_fricmatmult," time_fsample ",
267      & time_fsample
268 #endif
269       return
270       end  
271 c-------------------------------------------------------------------------------
272       subroutine velverlet_step(itime)
273 c-------------------------------------------------------------------------------
274 c  Perform a single velocity Verlet step; the time step can be rescaled if 
275 c  increments in accelerations exceed the threshold
276 c-------------------------------------------------------------------------------
277       implicit real*8 (a-h,o-z)
278       include 'DIMENSIONS'
279 #ifdef MPI
280       include 'mpif.h'
281       integer ierror,ierrcode
282 #endif
283       include 'COMMON.SETUP'
284       include 'COMMON.CONTROL'
285       include 'COMMON.VAR'
286       include 'COMMON.MD'
287 #ifndef LANG0
288       include 'COMMON.LANGEVIN'
289 #else
290       include 'COMMON.LANGEVIN.lang0'
291 #endif
292       include 'COMMON.CHAIN'
293       include 'COMMON.DERIV'
294       include 'COMMON.GEO'
295       include 'COMMON.LOCAL'
296       include 'COMMON.INTERACT'
297       include 'COMMON.IOUNITS'
298       include 'COMMON.NAMES'
299       include 'COMMON.TIME1'
300       include 'COMMON.MUCA'
301       double precision vcm(3),incr(3)
302       double precision cm(3),L(3)
303       integer ilen,count,rstcount
304       external ilen
305       character*50 tytul
306       integer maxcount_scale /20/
307       common /gucio/ cm
308       double precision stochforcvec(MAXRES6)
309       common /stochcalc/ stochforcvec
310       integer itime
311       logical scale
312       double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
313       double precision vtnp_(maxres6),vtnp_a(maxres6)
314 c
315       scale=.true.
316       icount_scale=0
317       if (lang.eq.1) then
318         call sddir_precalc
319       else if (lang.eq.2 .or. lang.eq.3) then
320 #ifndef LANG0
321         call stochastic_force(stochforcvec)
322 #else
323         write (iout,*) 
324      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
325 #ifdef MPI
326         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
327 #endif
328         stop
329 #endif
330       endif
331       itime_scal=0
332       do while (scale)
333         icount_scale=icount_scale+1
334         if (icount_scale.gt.maxcount_scale) then
335           write (iout,*) 
336      &      "ERROR: too many attempts at scaling down the time step. ",
337      &      "amax=",amax,"epdrift=",epdrift,
338      &      "damax=",damax,"edriftmax=",edriftmax,
339      &      "d_time=",d_time
340           call flush(iout)
341 #ifdef MPI
342           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
343 #endif
344           stop
345         endif
346 c First step of the velocity Verlet algorithm
347         if (lang.eq.2) then
348 #ifndef LANG0
349           call sd_verlet1
350 #endif
351         else if (lang.eq.3) then
352 #ifndef LANG0
353           call sd_verlet1_ciccotti
354 #endif
355         else if (lang.eq.1) then
356           call sddir_verlet1
357         else if (tnp1) then
358           call tnp1_step1
359         else if (tnp) then
360           call tnp_step1
361         else    
362           if (tnh) then
363             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
364             do i=0,2*nres
365              do j=1,3
366               d_t_old(j,i)=d_t_old(j,i)*scale_nh
367              enddo
368             enddo 
369           endif
370           call verlet1
371         endif     
372 c Build the chain from the newly calculated coordinates 
373         call chainbuild_cart
374         if (rattle) call rattle1
375         if (ntwe.ne.0) then
376         if (large.and. mod(itime,ntwe).eq.0) then
377           write (iout,*) "Cartesian and internal coordinates: step 1"
378           call cartprint
379           call intout
380           write (iout,*) "dC"
381           do i=0,nres
382             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
383      &      (dc(j,i+nres),j=1,3)
384           enddo
385           write (iout,*) "Accelerations"
386           do i=0,nres
387             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
388      &      (d_a(j,i+nres),j=1,3)
389           enddo
390           write (iout,*) "Velocities, step 1"
391           do i=0,nres
392             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
393      &      (d_t(j,i+nres),j=1,3)
394           enddo
395         endif
396         endif
397 #ifdef MPI
398         tt0 = MPI_Wtime()
399 #else
400         tt0 = tcpu()
401 #endif
402 c Calculate energy and forces
403         call zerograd
404         call etotal(potEcomp)
405 #ifdef TIMING_ENE
406 #ifdef MPI
407         t_etotal=t_etotal+MPI_Wtime()-tt0
408 #else
409         t_etotal=t_etotal+tcpu()-tt0
410 #endif
411 #endif
412         E_old=potE
413         potE=potEcomp(0)-potEcomp(20)
414         call cartgrad
415 c Get the new accelerations
416         call lagrangian
417 #ifdef MPI
418         t_enegrad=t_enegrad+MPI_Wtime()-tt0
419 #else
420         t_enegrad=t_enegrad+tcpu()-tt0
421 #endif
422 c Determine maximum acceleration and scale down the timestep if needed
423         call max_accel
424         amax=amax/(itime_scal+1)**2
425         call predict_edrift(epdrift)
426         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
427 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
428           scale=.true.
429           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
430      &      /dlog(2.0d0)+1
431           itime_scal=itime_scal+ifac_time
432 c          fac_time=dmin1(damax/amax,0.5d0)
433           fac_time=0.5d0**ifac_time
434           d_time=d_time*fac_time
435           if (lang.eq.2 .or. lang.eq.3) then 
436 #ifndef LANG0
437 c            write (iout,*) "Calling sd_verlet_setup: 1"
438 c Rescale the stochastic forces and recalculate or restore 
439 c the matrices of tinker integrator
440             if (itime_scal.gt.maxflag_stoch) then
441               if (large) write (iout,'(a,i5,a)') 
442      &         "Calculate matrices for stochastic step;",
443      &         " itime_scal ",itime_scal
444               if (lang.eq.2) then
445                 call sd_verlet_p_setup
446               else
447                 call sd_verlet_ciccotti_setup
448               endif
449               write (iout,'(2a,i3,a,i3,1h.)') 
450      &         "Warning: cannot store matrices for stochastic",
451      &         " integration because the index",itime_scal,
452      &         " is greater than",maxflag_stoch
453               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
454      &         " integration Langevin algorithm for better efficiency."
455             else if (flag_stoch(itime_scal)) then
456               if (large) write (iout,'(a,i5,a,l1)') 
457      &         "Restore matrices for stochastic step; itime_scal ",
458      &         itime_scal," flag ",flag_stoch(itime_scal)
459               do i=1,dimen3
460                 do j=1,dimen3
461                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
462                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
463                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
464                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
465                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
466                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
467                 enddo
468               enddo
469             else
470               if (large) write (iout,'(2a,i5,a,l1)') 
471      &         "Calculate & store matrices for stochastic step;",
472      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
473               if (lang.eq.2) then
474                 call sd_verlet_p_setup  
475               else
476                 call sd_verlet_ciccotti_setup
477               endif
478               flag_stoch(ifac_time)=.true.
479               do i=1,dimen3
480                 do j=1,dimen3
481                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
482                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
483                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
484                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
485                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
486                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
487                 enddo
488               enddo
489             endif
490             fac_time=1.0d0/dsqrt(fac_time)
491             do i=1,dimen3
492               stochforcvec(i)=fac_time*stochforcvec(i)
493             enddo
494 #endif
495           else if (lang.eq.1) then
496 c Rescale the accelerations due to stochastic forces
497             fac_time=1.0d0/dsqrt(fac_time)
498             do i=1,dimen3
499               d_as_work(i)=d_as_work(i)*fac_time
500             enddo
501           endif
502           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')  
503      &      "itime",itime," Timestep scaled down to ",
504      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
505         else 
506 c Second step of the velocity Verlet algorithm
507           if (lang.eq.2) then   
508 #ifndef LANG0
509             call sd_verlet2
510 #endif
511           else if (lang.eq.3) then
512 #ifndef LANG0
513             call sd_verlet2_ciccotti
514 #endif
515           else if (lang.eq.1) then
516             call sddir_verlet2
517           else if (tnp1) then
518             call tnp1_step2
519           else if (tnp) then
520             call tnp_step2
521           else
522             call verlet2
523             if (tnh) then
524               call kinetic(EK)
525               call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
526               do i=0,2*nres
527                do j=1,3
528                 d_t(j,i)=d_t(j,i)*scale_nh
529                enddo
530               enddo 
531             endif
532           endif                     
533           if (rattle) call rattle2
534           totT=totT+d_time
535           if (d_time.ne.d_time0) then
536             d_time=d_time0
537 #ifndef   LANG0
538             if (lang.eq.2 .or. lang.eq.3) then
539               if (large) write (iout,'(a)') 
540      &         "Restore original matrices for stochastic step"
541 c              write (iout,*) "Calling sd_verlet_setup: 2"
542 c Restore the matrices of tinker integrator if the time step has been restored
543               do i=1,dimen3
544                 do j=1,dimen3
545                   pfric_mat(i,j)=pfric0_mat(i,j,0)
546                   afric_mat(i,j)=afric0_mat(i,j,0)
547                   vfric_mat(i,j)=vfric0_mat(i,j,0)
548                   prand_mat(i,j)=prand0_mat(i,j,0)
549                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
550                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
551                 enddo
552               enddo
553             endif
554 #endif
555           endif
556           scale=.false.
557         endif
558       enddo
559 c Calculate the kinetic and the total energy and the kinetic temperature
560       if (tnp .or. tnp1) then 
561        do i=0,2*nres
562         do j=1,3
563           d_t_old(j,i)=d_t(j,i)
564           d_t(j,i)=d_t(j,i)/s_np
565         enddo
566        enddo 
567       endif
568       call kinetic(EK)
569       totE=EK+potE
570 c diagnostics
571 c      call kinetic1(EK1)
572 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
573 c end diagnostics
574 c Couple the system to Berendsen bath if needed
575       if (tbf .and. lang.eq.0) then
576         call verlet_bath
577       endif
578       kinetic_T=2.0d0/(dimen3*Rb)*EK
579 c Backup the coordinates, velocities, and accelerations
580       do i=0,2*nres
581         do j=1,3
582           dc_old(j,i)=dc(j,i)
583           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
584           d_a_old(j,i)=d_a(j,i)
585         enddo
586       enddo 
587       if (ntwe.ne.0) then
588       if (mod(itime,ntwe).eq.0) then
589
590        if(tnp .or. tnp1) then
591         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
592         H=(HNose1-H0)*s_np
593 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
594 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
595 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
596           hhh=h
597        endif
598
599        if(tnh) then
600         HNose1=Hnose_nh(EK,potE)
601         H=HNose1-H0
602         hhh=h
603 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
604        endif
605
606        if (large) then
607         itnp=0
608         do j=1,3
609          itnp=itnp+1
610          vtnp(itnp)=d_t(j,0)
611         enddo
612         do i=nnt,nct-1  
613          do j=1,3    
614           itnp=itnp+1
615           vtnp(itnp)=d_t(j,i)
616          enddo
617         enddo
618         do i=nnt,nct
619          if (itype(i).ne.10) then
620           inres=i+nres
621           do j=1,3  
622            itnp=itnp+1  
623            vtnp(itnp)=d_t(j,inres)
624           enddo
625          endif      
626         enddo 
627
628 c Transform velocities from UNRES coordinate space to cartesian and Gvec
629 c eigenvector space
630
631         do i=1,dimen3
632           vtnp_(i)=0.0d0
633           vtnp_a(i)=0.0d0
634           do j=1,dimen3
635             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
636             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
637           enddo
638           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
639         enddo
640
641         do i=1,dimen3
642          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
643         enddo
644
645         write (iout,*) "Velocities, step 2"
646         do i=0,nres
647           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
648      &      (d_t(j,i+nres),j=1,3)
649         enddo
650        endif
651       endif
652       endif
653       return
654       end
655 c-------------------------------------------------------------------------------
656       subroutine RESPA_step(itime)
657 c-------------------------------------------------------------------------------
658 c  Perform a single RESPA step.
659 c-------------------------------------------------------------------------------
660       implicit real*8 (a-h,o-z)
661       include 'DIMENSIONS'
662 #ifdef MPI
663       include 'mpif.h'
664       integer IERROR,ERRCODE
665 #endif
666       include 'COMMON.SETUP'
667       include 'COMMON.CONTROL'
668       include 'COMMON.VAR'
669       include 'COMMON.MD'
670 #ifndef LANG0
671       include 'COMMON.LANGEVIN'
672 #else
673       include 'COMMON.LANGEVIN.lang0'
674 #endif
675       include 'COMMON.CHAIN'
676       include 'COMMON.DERIV'
677       include 'COMMON.GEO'
678       include 'COMMON.LOCAL'
679       include 'COMMON.INTERACT'
680       include 'COMMON.IOUNITS'
681       include 'COMMON.NAMES'
682       include 'COMMON.TIME1'
683       double precision energia_short(0:n_ene),
684      & energia_long(0:n_ene)
685       double precision cm(3),L(3),vcm(3),incr(3)
686       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
687      & d_a_old0(3,0:maxres2)
688       integer ilen,count,rstcount
689       external ilen
690       character*50 tytul
691       integer maxcount_scale /10/
692       common /gucio/ cm,energia_short
693       double precision stochforcvec(MAXRES6)
694       common /stochcalc/ stochforcvec
695       integer itime
696       logical scale
697       double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
698       common /cipiszcze/ itt
699       itt=itime
700       if (ntwe.ne.0) then
701       if (large.and. mod(itime,ntwe).eq.0) then
702         write (iout,*) "***************** RESPA itime",itime
703         write (iout,*) "Cartesian and internal coordinates: step 0"
704 c        call cartprint
705         call pdbout(0.0d0,"cipiszcze",iout)
706         call intout
707         write (iout,*) "Accelerations from long-range forces"
708         do i=0,nres
709           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
710      &      (d_a(j,i+nres),j=1,3)
711         enddo
712         write (iout,*) "Velocities, step 0"
713         do i=0,nres
714           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
715      &      (d_t(j,i+nres),j=1,3)
716         enddo
717       endif
718       endif
719 c
720 c Perform the initial RESPA step (increment velocities)
721 c      write (iout,*) "*********************** RESPA ini"
722       if (tnp1) then
723           call tnp_respa_step1
724       else if (tnp) then
725           call tnp_respa_step1
726       else
727           if (tnh.and..not.xiresp) then
728             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
729             do i=0,2*nres
730              do j=1,3
731               d_t(j,i)=d_t(j,i)*scale_nh
732              enddo
733             enddo 
734           endif
735           call RESPA_vel
736       endif
737
738 cd       if(tnp .or. tnp1) then
739 cd        write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
740 cd       endif
741
742       if (ntwe.ne.0) then
743       if (mod(itime,ntwe).eq.0 .and. large) then
744         write (iout,*) "Velocities, end"
745         do i=0,nres
746           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
747      &      (d_t(j,i+nres),j=1,3)
748         enddo
749       endif
750       endif
751 c Compute the short-range forces
752 #ifdef MPI
753       tt0 =MPI_Wtime()
754 #else
755       tt0 = tcpu()
756 #endif
757 C 7/2/2009 commented out
758 c      call zerograd
759 c      call etotal_short(energia_short)
760       if (tnp.or.tnp1) potE=energia_short(0)
761 c      call cartgrad
762 c      call lagrangian
763 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
764         do i=0,2*nres
765           do j=1,3
766             d_a(j,i)=d_a_short(j,i)
767           enddo
768         enddo
769       if (ntwe.ne.0) then
770       if (large.and. mod(itime,ntwe).eq.0) then
771         write (iout,*) "energia_short",energia_short(0)
772         write (iout,*) "Accelerations from short-range forces"
773         do i=0,nres
774           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
775      &      (d_a(j,i+nres),j=1,3)
776         enddo
777       endif
778       endif
779 #ifdef MPI
780         t_enegrad=t_enegrad+MPI_Wtime()-tt0
781 #else
782         t_enegrad=t_enegrad+tcpu()-tt0
783 #endif
784       do i=0,2*nres
785         do j=1,3
786           dc_old(j,i)=dc(j,i)
787           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
788           d_a_old(j,i)=d_a(j,i)
789         enddo
790       enddo 
791 c 6/30/08 A-MTS: attempt at increasing the split number
792       do i=0,2*nres
793         do j=1,3
794           dc_old0(j,i)=dc_old(j,i)
795           d_t_old0(j,i)=d_t_old(j,i)
796           d_a_old0(j,i)=d_a_old(j,i)
797         enddo
798       enddo 
799       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
800       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
801 c
802       scale=.true.
803       d_time0=d_time
804       do while (scale)
805
806       scale=.false.
807 c      write (iout,*) "itime",itime," ntime_split",ntime_split
808 c Split the time step
809       d_time=d_time0/ntime_split 
810 c Perform the short-range RESPA steps (velocity Verlet increments of
811 c positions and velocities using short-range forces)
812 c      write (iout,*) "*********************** RESPA split"
813       do itsplit=1,ntime_split
814         if (lang.eq.1) then
815           call sddir_precalc
816         else if (lang.eq.2 .or. lang.eq.3) then
817 #ifndef LANG0
818           call stochastic_force(stochforcvec)
819 #else
820           write (iout,*) 
821      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
822 #ifdef MPI
823           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
824 #endif
825           stop
826 #endif
827         endif
828 c First step of the velocity Verlet algorithm
829         if (lang.eq.2) then
830 #ifndef LANG0
831           call sd_verlet1
832 #endif
833         else if (lang.eq.3) then
834 #ifndef LANG0
835           call sd_verlet1_ciccotti
836 #endif
837         else if (lang.eq.1) then
838           call sddir_verlet1
839         else if (tnp1) then
840           call tnp1_respa_i_step1
841         else if (tnp) then
842           call tnp_respa_i_step1
843         else
844           if (tnh.and.xiresp) then
845             call kinetic(EK)
846             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
847             do i=0,2*nres
848              do j=1,3
849               d_t_old(j,i)=d_t_old(j,i)*scale_nh
850              enddo
851             enddo 
852 cd            write(iout,*) "SSS1",itsplit,EK,scale_nh
853           endif
854           call verlet1
855         endif
856 c Build the chain from the newly calculated coordinates 
857         call chainbuild_cart
858         if (rattle) call rattle1
859         if (ntwe.ne.0) then
860         if (large.and. mod(itime,ntwe).eq.0) then
861           write (iout,*) "***** ITSPLIT",itsplit
862           write (iout,*) "Cartesian and internal coordinates: step 1"
863           call pdbout(0.0d0,"cipiszcze",iout)
864 c          call cartprint
865           call intout
866           write (iout,*) "Velocities, step 1"
867           do i=0,nres
868             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
869      &        (d_t(j,i+nres),j=1,3)
870           enddo
871         endif
872         endif
873 #ifdef MPI
874         tt0 = MPI_Wtime()
875 #else
876         tt0 = tcpu()
877 #endif
878 c Calculate energy and forces
879         call zerograd
880         call etotal_short(energia_short)
881         E_old=potE
882         potE=energia_short(0)
883 #ifdef TIMING_ENE
884 #ifdef MPI
885         t_eshort=t_eshort+MPI_Wtime()-tt0
886 #else
887         t_eshort=t_eshort+tcpu()-tt0
888 #endif
889 #endif
890         call cartgrad
891 c Get the new accelerations
892         call lagrangian
893 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
894         do i=0,2*nres
895           do j=1,3
896             d_a_short(j,i)=d_a(j,i)
897           enddo
898         enddo
899         if (ntwe.ne.0) then
900         if (large.and. mod(itime,ntwe).eq.0) then
901           write (iout,*)"energia_short",energia_short(0)
902           write (iout,*) "Accelerations from short-range forces"
903           do i=0,nres
904             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
905      &        (d_a(j,i+nres),j=1,3)
906           enddo
907         endif
908         endif
909 c 6/30/08 A-MTS
910 c Determine maximum acceleration and scale down the timestep if needed
911         call max_accel
912         amax=amax/ntime_split**2
913         call predict_edrift(epdrift)
914         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) 
915      &   write (iout,*) "amax",amax," damax",damax,
916      &   " epdrift",epdrift," epdriftmax",epdriftmax
917 c Exit loop and try with increased split number if the change of
918 c acceleration is too big
919         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
920           if (ntime_split.lt.maxtime_split) then
921             scale=.true.
922             ntime_split=ntime_split*2
923             do i=0,2*nres
924               do j=1,3
925                 dc_old(j,i)=dc_old0(j,i)
926                 d_t_old(j,i)=d_t_old0(j,i)
927                 d_a_old(j,i)=d_a_old0(j,i)
928               enddo
929             enddo 
930             write (iout,*) "acceleration/energy drift too large",amax,
931      &      epdrift," split increased to ",ntime_split," itime",itime,
932      &       " itsplit",itsplit
933             exit
934           else
935             write (iout,*) 
936      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
937      &       maxtime_split,
938      &      " already reached!!! Trying to carry on!"
939           endif
940         endif
941 #ifdef MPI
942         t_enegrad=t_enegrad+MPI_Wtime()-tt0
943 #else
944         t_enegrad=t_enegrad+tcpu()-tt0
945 #endif
946 c Second step of the velocity Verlet algorithm
947         if (lang.eq.2) then
948 #ifndef LANG0
949           call sd_verlet2
950 #endif
951         else if (lang.eq.3) then
952 #ifndef LANG0
953           call sd_verlet2_ciccotti
954 #endif
955         else if (lang.eq.1) then
956           call sddir_verlet2
957         else if (tnp1) then
958             call tnp1_respa_i_step2
959         else if (tnp) then
960             call tnp_respa_i_step2
961         else
962           call verlet2
963           if (tnh.and.xiresp) then
964             call kinetic(EK)
965             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
966             do i=0,2*nres
967              do j=1,3
968               d_t(j,i)=d_t(j,i)*scale_nh
969              enddo
970             enddo 
971 cd            write(iout,*) "SSS2",itsplit,EK,scale_nh
972           endif
973         endif
974         if (rattle) call rattle2
975 c Backup the coordinates, velocities, and accelerations
976         if (tnp .or. tnp1) then 
977          do i=0,2*nres
978           do j=1,3
979             d_t_old(j,i)=d_t(j,i)
980             if (tnp) d_t(j,i)=d_t(j,i)/s_np
981             if (tnp1) d_t(j,i)=d_t(j,i)/s_np
982           enddo
983          enddo 
984         endif
985
986         do i=0,2*nres
987           do j=1,3
988             dc_old(j,i)=dc(j,i)
989             if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
990             d_a_old(j,i)=d_a(j,i)
991           enddo
992         enddo 
993       enddo
994
995       enddo ! while scale
996
997 c Restore the time step
998       d_time=d_time0
999 c Compute long-range forces
1000 #ifdef MPI
1001       tt0 =MPI_Wtime()
1002 #else
1003       tt0 = tcpu()
1004 #endif
1005       call zerograd
1006       call etotal_long(energia_long)
1007       E_long=energia_long(0)
1008       potE=energia_short(0)+energia_long(0)
1009 #ifdef TIMING_ENE
1010 #ifdef MPI
1011         t_elong=t_elong+MPI_Wtime()-tt0
1012 #else
1013         t_elong=t_elong+tcpu()-tt0
1014 #endif
1015 #endif
1016       call cartgrad
1017       call lagrangian
1018 #ifdef MPI
1019         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1020 #else
1021         t_enegrad=t_enegrad+tcpu()-tt0
1022 #endif
1023 c Compute accelerations from long-range forces
1024       if (ntwe.ne.0) then
1025       if (large.and. mod(itime,ntwe).eq.0) then
1026         write (iout,*) "energia_long",energia_long(0)
1027         write (iout,*) "Cartesian and internal coordinates: step 2"
1028 c        call cartprint
1029         call pdbout(0.0d0,"cipiszcze",iout)
1030         call intout
1031         write (iout,*) "Accelerations from long-range forces"
1032         do i=0,nres
1033           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1034      &      (d_a(j,i+nres),j=1,3)
1035         enddo
1036         write (iout,*) "Velocities, step 2"
1037         do i=0,nres
1038           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1039      &      (d_t(j,i+nres),j=1,3)
1040         enddo
1041       endif
1042       endif
1043 c Compute the final RESPA step (increment velocities)
1044 c      write (iout,*) "*********************** RESPA fin"
1045       if (tnp1) then
1046           call tnp_respa_step2
1047       else if (tnp) then
1048           call tnp_respa_step2
1049       else
1050           call RESPA_vel
1051           if (tnh.and..not.xiresp) then
1052             call kinetic(EK)
1053             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1054             do i=0,2*nres
1055              do j=1,3
1056               d_t(j,i)=d_t(j,i)*scale_nh
1057              enddo
1058             enddo 
1059           endif
1060       endif
1061
1062         if (tnp .or. tnp1) then 
1063          do i=0,2*nres
1064           do j=1,3
1065             d_t(j,i)=d_t_old(j,i)/s_np
1066           enddo
1067          enddo 
1068         endif
1069
1070 c Compute the complete potential energy
1071       do i=0,n_ene
1072         potEcomp(i)=energia_short(i)+energia_long(i)
1073       enddo
1074       potE=potEcomp(0)-potEcomp(20)
1075 c      potE=energia_short(0)+energia_long(0)
1076       totT=totT+d_time
1077 c Calculate the kinetic and the total energy and the kinetic temperature
1078       call kinetic(EK)
1079       totE=EK+potE
1080 c Couple the system to Berendsen bath if needed
1081       if (tbf .and. lang.eq.0) then
1082         call verlet_bath
1083       endif
1084       kinetic_T=2.0d0/(dimen3*Rb)*EK
1085 c Backup the coordinates, velocities, and accelerations
1086       if (ntwe.ne.0) then
1087       if (mod(itime,ntwe).eq.0 .and. large) then
1088         write (iout,*) "Velocities, end"
1089         do i=0,nres
1090           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1091      &      (d_t(j,i+nres),j=1,3)
1092         enddo
1093       endif
1094
1095       if (mod(itime,ntwe).eq.0) then
1096
1097        if(tnp .or. tnp1) then
1098 #ifndef G77
1099         write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
1100      &                          E_long,energia_short(0)
1101 #else
1102         write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1103      &                          E_long,energia_short(0)
1104 #endif
1105         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1106         H=(HNose1-H0)*s_np
1107 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1108 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
1109 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1110           hhh=h
1111 cd        write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1112        endif
1113
1114        if(tnh) then
1115         HNose1=Hnose_nh(EK,potE)
1116         H=HNose1-H0
1117 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1118         hhh=h
1119        endif
1120
1121
1122        if (large) then
1123        itnp=0
1124        do j=1,3
1125         itnp=itnp+1
1126         vtnp(itnp)=d_t(j,0)
1127        enddo
1128        do i=nnt,nct-1   
1129         do j=1,3    
1130           itnp=itnp+1
1131           vtnp(itnp)=d_t(j,i)
1132         enddo
1133        enddo
1134        do i=nnt,nct
1135         if (itype(i).ne.10) then
1136           inres=i+nres
1137           do j=1,3  
1138            itnp=itnp+1  
1139            vtnp(itnp)=d_t(j,inres)
1140           enddo
1141         endif      
1142        enddo 
1143
1144 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1145 c eigenvector space
1146
1147         do i=1,dimen3
1148           vtnp_(i)=0.0d0
1149           vtnp_a(i)=0.0d0
1150           do j=1,dimen3
1151             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1152             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1153           enddo
1154           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1155         enddo
1156
1157         do i=1,dimen3
1158          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1159         enddo
1160
1161        endif
1162       endif
1163       endif
1164
1165
1166       return
1167       end
1168 c---------------------------------------------------------------------
1169       subroutine RESPA_vel
1170 c  First and last RESPA step (incrementing velocities using long-range
1171 c  forces).
1172       implicit real*8 (a-h,o-z)
1173       include 'DIMENSIONS'
1174       include 'COMMON.CONTROL'
1175       include 'COMMON.VAR'
1176       include 'COMMON.MD'
1177       include 'COMMON.CHAIN'
1178       include 'COMMON.DERIV'
1179       include 'COMMON.GEO'
1180       include 'COMMON.LOCAL'
1181       include 'COMMON.INTERACT'
1182       include 'COMMON.IOUNITS'
1183       include 'COMMON.NAMES'
1184       do j=1,3
1185         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1186       enddo
1187       do i=nnt,nct-1
1188         do j=1,3
1189           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1190         enddo
1191       enddo
1192       do i=nnt,nct
1193         if (itype(i).ne.10) then
1194           inres=i+nres
1195           do j=1,3
1196             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1197           enddo
1198         endif
1199       enddo 
1200       return
1201       end
1202 c-----------------------------------------------------------------
1203       subroutine verlet1
1204 c Applying velocity Verlet algorithm - step 1 to coordinates
1205       implicit real*8 (a-h,o-z)
1206       include 'DIMENSIONS'
1207       include 'COMMON.CONTROL'
1208       include 'COMMON.VAR'
1209       include 'COMMON.MD'
1210       include 'COMMON.CHAIN'
1211       include 'COMMON.DERIV'
1212       include 'COMMON.GEO'
1213       include 'COMMON.LOCAL'
1214       include 'COMMON.INTERACT'
1215       include 'COMMON.IOUNITS'
1216       include 'COMMON.NAMES'
1217       double precision adt,adt2
1218         
1219 #ifdef DEBUG
1220       write (iout,*) "VELVERLET1 START: DC"
1221       do i=0,nres
1222         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1223      &   (dc(j,i+nres),j=1,3)
1224       enddo 
1225 #endif
1226       do j=1,3
1227         adt=d_a_old(j,0)*d_time
1228         adt2=0.5d0*adt
1229         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1230         d_t_new(j,0)=d_t_old(j,0)+adt2
1231         d_t(j,0)=d_t_old(j,0)+adt
1232       enddo
1233       do i=nnt,nct-1    
1234         do j=1,3    
1235           adt=d_a_old(j,i)*d_time
1236           adt2=0.5d0*adt
1237           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1238           d_t_new(j,i)=d_t_old(j,i)+adt2
1239           d_t(j,i)=d_t_old(j,i)+adt
1240         enddo
1241       enddo
1242       do i=nnt,nct
1243         if (itype(i).ne.10) then
1244           inres=i+nres
1245           do j=1,3    
1246             adt=d_a_old(j,inres)*d_time
1247             adt2=0.5d0*adt
1248             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1249             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1250             d_t(j,inres)=d_t_old(j,inres)+adt
1251           enddo
1252         endif      
1253       enddo 
1254 #ifdef DEBUG
1255       write (iout,*) "VELVERLET1 END: DC"
1256       do i=0,nres
1257         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1258      &   (dc(j,i+nres),j=1,3)
1259       enddo 
1260 #endif
1261       return
1262       end
1263 c---------------------------------------------------------------------
1264       subroutine verlet2
1265 c  Step 2 of the velocity Verlet algorithm: update velocities
1266       implicit real*8 (a-h,o-z)
1267       include 'DIMENSIONS'
1268       include 'COMMON.CONTROL'
1269       include 'COMMON.VAR'
1270       include 'COMMON.MD'
1271       include 'COMMON.CHAIN'
1272       include 'COMMON.DERIV'
1273       include 'COMMON.GEO'
1274       include 'COMMON.LOCAL'
1275       include 'COMMON.INTERACT'
1276       include 'COMMON.IOUNITS'
1277       include 'COMMON.NAMES'
1278       do j=1,3
1279         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1280       enddo
1281       do i=nnt,nct-1
1282         do j=1,3
1283           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1284         enddo
1285       enddo
1286       do i=nnt,nct
1287         if (itype(i).ne.10) then
1288           inres=i+nres
1289           do j=1,3
1290             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1291           enddo
1292         endif
1293       enddo 
1294       return
1295       end
1296 c-----------------------------------------------------------------
1297       subroutine sddir_precalc
1298 c Applying velocity Verlet algorithm - step 1 to coordinates        
1299       implicit real*8 (a-h,o-z)
1300       include 'DIMENSIONS'
1301 #ifdef MPI
1302       include 'mpif.h'
1303 #endif
1304       include 'COMMON.CONTROL'
1305       include 'COMMON.VAR'
1306       include 'COMMON.MD'
1307 #ifndef LANG0
1308       include 'COMMON.LANGEVIN'
1309 #else
1310       include 'COMMON.LANGEVIN.lang0'
1311 #endif
1312       include 'COMMON.CHAIN'
1313       include 'COMMON.DERIV'
1314       include 'COMMON.GEO'
1315       include 'COMMON.LOCAL'
1316       include 'COMMON.INTERACT'
1317       include 'COMMON.IOUNITS'
1318       include 'COMMON.NAMES'
1319       include 'COMMON.TIME1'
1320       double precision stochforcvec(MAXRES6)
1321       common /stochcalc/ stochforcvec
1322 c
1323 c Compute friction and stochastic forces
1324 c
1325 #ifdef MPI
1326       time00=MPI_Wtime()
1327 #else
1328       time00=tcpu()
1329 #endif
1330       call friction_force
1331 #ifdef MPI
1332       time_fric=time_fric+MPI_Wtime()-time00
1333       time00=MPI_Wtime()
1334 #else
1335       time_fric=time_fric+tcpu()-time00
1336       time00=tcpu()
1337 #endif
1338       call stochastic_force(stochforcvec) 
1339 #ifdef MPI
1340       time_stoch=time_stoch+MPI_Wtime()-time00
1341 #else
1342       time_stoch=time_stoch+tcpu()-time00
1343 #endif
1344 c
1345 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1346 c forces (d_as_work)
1347 c
1348       call ginv_mult(fric_work, d_af_work)
1349       call ginv_mult(stochforcvec, d_as_work)
1350       return
1351       end
1352 c---------------------------------------------------------------------
1353       subroutine sddir_verlet1
1354 c Applying velocity Verlet algorithm - step 1 to velocities        
1355       implicit real*8 (a-h,o-z)
1356       include 'DIMENSIONS'
1357       include 'COMMON.CONTROL'
1358       include 'COMMON.VAR'
1359       include 'COMMON.MD'
1360 #ifndef LANG0
1361       include 'COMMON.LANGEVIN'
1362 #else
1363       include 'COMMON.LANGEVIN.lang0'
1364 #endif
1365       include 'COMMON.CHAIN'
1366       include 'COMMON.DERIV'
1367       include 'COMMON.GEO'
1368       include 'COMMON.LOCAL'
1369       include 'COMMON.INTERACT'
1370       include 'COMMON.IOUNITS'
1371       include 'COMMON.NAMES'
1372 c Revised 3/31/05 AL: correlation between random contributions to 
1373 c position and velocity increments included.
1374       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1375       double precision adt,adt2
1376 c
1377 c Add the contribution from BOTH friction and stochastic force to the
1378 c coordinates, but ONLY the contribution from the friction forces to velocities
1379 c
1380       do j=1,3
1381         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1382         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1383         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1384         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1385         d_t(j,0)=d_t_old(j,0)+adt
1386       enddo
1387       ind=3
1388       do i=nnt,nct-1    
1389         do j=1,3    
1390           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1391           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1392           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1393           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1394           d_t(j,i)=d_t_old(j,i)+adt
1395         enddo
1396         ind=ind+3
1397       enddo
1398       do i=nnt,nct
1399         if (itype(i).ne.10) then
1400           inres=i+nres
1401           do j=1,3    
1402             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1403             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1404             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1405             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1406             d_t(j,inres)=d_t_old(j,inres)+adt
1407           enddo
1408           ind=ind+3
1409         endif      
1410       enddo 
1411       return
1412       end
1413 c---------------------------------------------------------------------
1414       subroutine sddir_verlet2
1415 c  Calculating the adjusted velocities for accelerations
1416       implicit real*8 (a-h,o-z)
1417       include 'DIMENSIONS'
1418       include 'COMMON.CONTROL'
1419       include 'COMMON.VAR'
1420       include 'COMMON.MD'
1421 #ifndef LANG0
1422       include 'COMMON.LANGEVIN'
1423 #else
1424       include 'COMMON.LANGEVIN.lang0'
1425 #endif
1426       include 'COMMON.CHAIN'
1427       include 'COMMON.DERIV'
1428       include 'COMMON.GEO'
1429       include 'COMMON.LOCAL'
1430       include 'COMMON.INTERACT'
1431       include 'COMMON.IOUNITS'
1432       include 'COMMON.NAMES'
1433       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1434       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1435 c Revised 3/31/05 AL: correlation between random contributions to 
1436 c position and velocity increments included.
1437 c The correlation coefficients are calculated at low-friction limit.
1438 c Also, friction forces are now not calculated with new velocities.
1439
1440 c      call friction_force
1441       call stochastic_force(stochforcvec) 
1442 c
1443 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1444 c forces (d_as_work)
1445 c
1446       call ginv_mult(stochforcvec, d_as_work1)
1447
1448 c
1449 c Update velocities
1450 c
1451       do j=1,3
1452         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1453      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1454       enddo
1455       ind=3
1456       do i=nnt,nct-1
1457         do j=1,3
1458           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1459      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1460         enddo
1461         ind=ind+3
1462       enddo
1463       do i=nnt,nct
1464         if (itype(i).ne.10) then
1465           inres=i+nres
1466           do j=1,3
1467             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1468      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1469      &       +cos60*d_as_work1(ind+j))*d_time
1470           enddo
1471           ind=ind+3
1472         endif
1473       enddo 
1474       return
1475       end
1476 c---------------------------------------------------------------------
1477       subroutine max_accel
1478 c
1479 c Find the maximum difference in the accelerations of the the sites
1480 c at the beginning and the end of the time step.
1481 c
1482       implicit real*8 (a-h,o-z)
1483       include 'DIMENSIONS'
1484       include 'COMMON.CONTROL'
1485       include 'COMMON.VAR'
1486       include 'COMMON.MD'
1487       include 'COMMON.CHAIN'
1488       include 'COMMON.DERIV'
1489       include 'COMMON.GEO'
1490       include 'COMMON.LOCAL'
1491       include 'COMMON.INTERACT'
1492       include 'COMMON.IOUNITS'
1493       double precision aux(3),accel(3),accel_old(3),dacc
1494       do j=1,3
1495 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1496          accel_old(j)=d_a_old(j,0)
1497          accel(j)=d_a(j,0)
1498       enddo 
1499       amax=0.0d0
1500       do i=nnt,nct
1501 c Backbone
1502         if (i.lt.nct) then
1503 c 7/3/08 changed to asymmetric difference
1504           do j=1,3
1505 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1506             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1507             accel(j)=accel(j)+0.5d0*d_a(j,i)
1508 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1509             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1510               dacc=dabs(accel(j)-accel_old(j))
1511               if (dacc.gt.amax) amax=dacc
1512             endif
1513           enddo
1514         endif
1515       enddo
1516 c Side chains
1517       do j=1,3
1518 c        accel(j)=aux(j)
1519         accel_old(j)=d_a_old(j,0)
1520         accel(j)=d_a(j,0)
1521       enddo
1522       if (nnt.eq.2) then
1523         do j=1,3
1524           accel_old(j)=accel_old(j)+d_a_old(j,1)
1525           accel(j)=accel(j)+d_a(j,1)
1526         enddo
1527       endif
1528       do i=nnt,nct
1529         if (itype(i).ne.10) then
1530           do j=1,3 
1531 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1532             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1533             accel(j)=accel(j)+d_a(j,i+nres)
1534           enddo
1535         endif
1536         do j=1,3
1537 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1538           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1539             dacc=dabs(accel(j)-accel_old(j))
1540             if (dacc.gt.amax) amax=dacc
1541           endif
1542         enddo
1543         do j=1,3
1544           accel_old(j)=accel_old(j)+d_a_old(j,i)
1545           accel(j)=accel(j)+d_a(j,i)
1546 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1547         enddo
1548       enddo
1549       return
1550       end       
1551 c---------------------------------------------------------------------
1552       subroutine predict_edrift(epdrift)
1553 c
1554 c Predict the drift of the potential energy
1555 c
1556       implicit real*8 (a-h,o-z)
1557       include 'DIMENSIONS'
1558       include 'COMMON.CONTROL'
1559       include 'COMMON.VAR'
1560       include 'COMMON.MD'
1561       include 'COMMON.CHAIN'
1562       include 'COMMON.DERIV'
1563       include 'COMMON.GEO'
1564       include 'COMMON.LOCAL'
1565       include 'COMMON.INTERACT'
1566       include 'COMMON.IOUNITS'
1567       include 'COMMON.MUCA'
1568       double precision epdrift,epdriftij
1569 c Drift of the potential energy
1570       epdrift=0.0d0
1571       do i=nnt,nct
1572 c Backbone
1573         if (i.lt.nct) then
1574           do j=1,3
1575             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1576             if (lmuca) epdriftij=epdriftij*factor
1577 c            write (iout,*) "back",i,j,epdriftij
1578             if (epdriftij.gt.epdrift) epdrift=epdriftij 
1579           enddo
1580         endif
1581 c Side chains
1582         if (itype(i).ne.10) then
1583           do j=1,3 
1584             epdriftij=
1585      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1586             if (lmuca) epdriftij=epdriftij*factor
1587 c            write (iout,*) "side",i,j,epdriftij
1588             if (epdriftij.gt.epdrift) epdrift=epdriftij
1589           enddo
1590         endif
1591       enddo
1592       epdrift=0.5d0*epdrift*d_time*d_time
1593 c      write (iout,*) "epdrift",epdrift
1594       return
1595       end       
1596 c-----------------------------------------------------------------------
1597       subroutine verlet_bath
1598 c
1599 c  Coupling to the thermostat by using the Berendsen algorithm
1600 c
1601       implicit real*8 (a-h,o-z)
1602       include 'DIMENSIONS'
1603       include 'COMMON.CONTROL'
1604       include 'COMMON.VAR'
1605       include 'COMMON.MD'
1606       include 'COMMON.CHAIN'
1607       include 'COMMON.DERIV'
1608       include 'COMMON.GEO'
1609       include 'COMMON.LOCAL'
1610       include 'COMMON.INTERACT'
1611       include 'COMMON.IOUNITS'
1612       include 'COMMON.NAMES'
1613       double precision T_half,fact
1614
1615       T_half=2.0d0/(dimen3*Rb)*EK
1616       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1617 c      write(iout,*) "T_half", T_half
1618 c      write(iout,*) "EK", EK
1619 c      write(iout,*) "fact", fact                               
1620       do j=1,3
1621         d_t(j,0)=fact*d_t(j,0)
1622       enddo
1623       do i=nnt,nct-1
1624         do j=1,3
1625           d_t(j,i)=fact*d_t(j,i)
1626         enddo
1627       enddo
1628       do i=nnt,nct
1629         if (itype(i).ne.10) then
1630           inres=i+nres
1631           do j=1,3
1632             d_t(j,inres)=fact*d_t(j,inres)
1633           enddo
1634         endif
1635       enddo 
1636       return
1637       end
1638 c---------------------------------------------------------
1639       subroutine init_MD
1640 c  Set up the initial conditions of a MD simulation
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643 #ifdef MP
1644       include 'mpif.h'
1645       character*16 form
1646       integer IERROR,ERRCODE
1647 #endif
1648       include 'COMMON.SETUP'
1649       include 'COMMON.CONTROL'
1650       include 'COMMON.VAR'
1651       include 'COMMON.MD'
1652 #ifndef LANG0
1653       include 'COMMON.LANGEVIN'
1654 #else
1655       include 'COMMON.LANGEVIN.lang0'
1656 #endif
1657       include 'COMMON.CHAIN'
1658       include 'COMMON.DERIV'
1659       include 'COMMON.GEO'
1660       include 'COMMON.LOCAL'
1661       include 'COMMON.INTERACT'
1662       include 'COMMON.IOUNITS'
1663       include 'COMMON.NAMES'
1664       include 'COMMON.REMD'
1665       real*8 energia_long(0:n_ene),
1666      &  energia_short(0:n_ene),vcm(3),incr(3),E_short
1667       double precision cm(3),L(3),xv,sigv,lowb,highb
1668       double precision varia(maxvar)
1669       character*256 qstr
1670       integer ilen
1671       external ilen
1672       character*50 tytul
1673       logical file_exist
1674       common /gucio/ cm
1675       d_time0=d_time
1676 c      write(iout,*) "d_time", d_time
1677 c Compute the standard deviations of stochastic forces for Langevin dynamics
1678 c if the friction coefficients do not depend on surface area
1679       if (lang.gt.0 .and. .not.surfarea) then
1680         do i=nnt,nct-1
1681           stdforcp(i)=stdfp*dsqrt(gamp)
1682         enddo
1683         do i=nnt,nct
1684           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1685         enddo
1686       endif
1687 c Open the pdb file for snapshotshots
1688 #ifdef MPI
1689       if(mdpdb) then
1690         if (ilen(tmpdir).gt.0) 
1691      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1692      &      liczba(:ilen(liczba))//".pdb")
1693         open(ipdb,
1694      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1695      &  //".pdb")
1696       else
1697 #ifdef NOXDR
1698         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
1699      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1700      &      liczba(:ilen(liczba))//".x")
1701         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1702      &  //".x"
1703 #else
1704         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
1705      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1706      &      liczba(:ilen(liczba))//".cx")
1707         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1708      &  //".cx"
1709 #endif
1710       endif
1711 #else
1712       if(mdpdb) then
1713          if (ilen(tmpdir).gt.0) 
1714      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1715          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1716       else
1717          if (ilen(tmpdir).gt.0) 
1718      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1719          cartname=prefix(:ilen(prefix))//"_MD.cx"
1720       endif
1721 #endif
1722       if (usampl) then
1723         write (qstr,'(256(1h ))')
1724         ipos=1
1725         do i=1,nfrag
1726           iq = qinfrag(i,iset)*10
1727           iw = wfrag(i,iset)/100
1728           if (iw.gt.0) then
1729             if(me.eq.king.or..not.out1file)
1730      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1731             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1732             ipos=ipos+7
1733           endif
1734         enddo
1735         do i=1,npair
1736           iq = qinpair(i,iset)*10
1737           iw = wpair(i,iset)/100
1738           if (iw.gt.0) then
1739             if(me.eq.king.or..not.out1file)
1740      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1741             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1742             ipos=ipos+7
1743           endif
1744         enddo
1745 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1746 #ifdef NOXDR
1747 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1748 #else
1749 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1750 #endif
1751 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1752       endif
1753       icg=1
1754       if (rest) then
1755        if (restart1file) then
1756          if (me.eq.king)
1757      &     inquire(file=mremd_rst_name,exist=file_exist)
1758            write (*,*) me," Before broadcast: file_exist",file_exist
1759 #ifdef MPI
1760          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1761      &          IERR)
1762          write (*,*) me," After broadcast: file_exist",file_exist
1763 #endif
1764 c        inquire(file=mremd_rst_name,exist=file_exist)
1765         if(me.eq.king.or..not.out1file)
1766      &   write(iout,*) "Initial state read by master and distributed"
1767        else
1768          if (ilen(tmpdir).gt.0)
1769      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1770      &      //liczba(:ilen(liczba))//'.rst')
1771         inquire(file=rest2name,exist=file_exist)
1772        endif
1773        if(file_exist) then
1774          if(.not.restart1file) then
1775            if(me.eq.king.or..not.out1file)
1776      &      write(iout,*) "Initial state will be read from file ",
1777      &      rest2name(:ilen(rest2name))
1778            call readrst
1779          endif  
1780          call rescale_weights(t_bath)
1781        else
1782         if(me.eq.king.or..not.out1file)then
1783          if (restart1file) then
1784           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1785      &       " does not exist"
1786          else
1787           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1788      &       " does not exist"
1789          endif
1790          write(iout,*) "Initial velocities randomly generated"
1791         endif
1792         call random_vel
1793         totT=0.0d0
1794        endif
1795       else
1796 c Generate initial velocities
1797         if(me.eq.king.or..not.out1file)
1798      &   write(iout,*) "Initial velocities randomly generated"
1799         call random_vel
1800         totT=0.0d0
1801       endif
1802 c      rest2name = prefix(:ilen(prefix))//'.rst'
1803       if(me.eq.king.or..not.out1file)then
1804        write (iout,*) "Initial velocities"
1805        do i=0,nres
1806          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1807      &   (d_t(j,i+nres),j=1,3)
1808        enddo                     
1809        call flush(iout)
1810 c  Zeroing the total angular momentum of the system
1811        write(iout,*) "Calling the zero-angular 
1812      & momentum subroutine"
1813       endif
1814       call inertia_tensor  
1815 c  Getting the potential energy and forces and velocities and accelerations
1816       call vcm_vel(vcm)
1817 c      write (iout,*) "velocity of the center of the mass:"
1818 c      write (iout,*) (vcm(j),j=1,3)
1819       do j=1,3
1820         d_t(j,0)=d_t(j,0)-vcm(j)
1821       enddo
1822 c Removing the velocity of the center of mass
1823       call vcm_vel(vcm)
1824       if(me.eq.king.or..not.out1file)then
1825        write (iout,*) "vcm right after adjustment:"
1826        write (iout,*) (vcm(j),j=1,3) 
1827        call flush(iout)
1828       endif
1829       if (.not.rest) then               
1830          call chainbuild
1831          if(iranconf.ne.0) then
1832           if (overlapsc) then 
1833            print *, 'Calling OVERLAP_SC'
1834            call overlap_sc(fail)
1835           endif 
1836
1837           if (searchsc) then 
1838            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1839            print *,'SC_move',nft_sc,etot
1840            if(me.eq.king.or..not.out1file)
1841      &      write(iout,*) 'SC_move',nft_sc,etot
1842           endif 
1843
1844           if(dccart)then
1845            print *, 'Calling MINIM_DC'
1846            call minim_dc(etot,iretcode,nfun)
1847           else
1848            call geom_to_var(nvar,varia)
1849            print *,'Calling MINIMIZE.'
1850            call minimize(etot,varia,iretcode,nfun)
1851            call var_to_geom(nvar,varia)
1852           endif
1853           if(me.eq.king.or..not.out1file)
1854      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1855          endif
1856       endif       
1857       call chainbuild_cart
1858       call kinetic(EK)
1859       if (tbf) then
1860         call verlet_bath(EK)
1861       endif       
1862       kinetic_T=2.0d0/(dimen3*Rb)*EK
1863       if(me.eq.king.or..not.out1file)then
1864        call cartprint
1865        call intout
1866       endif
1867 #ifdef MPI
1868       tt0=MPI_Wtime()
1869 #else
1870       tt0=tcpu()
1871 #endif
1872       call zerograd
1873       call etotal(potEcomp)
1874 #ifdef TIMING_ENE
1875 #ifdef MPI
1876       t_etotal=t_etotal+MPI_Wtime()-tt0
1877 #else
1878       t_etotal=t_etotal+tcpu()-tt0
1879 #endif
1880 #endif
1881       potE=potEcomp(0)
1882
1883       if(tnp .or. tnp1) then
1884        s_np=1.0
1885        pi_np=0.0
1886        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1887        H0=Hnose1
1888        write(iout,*) 'H0= ',H0
1889       endif
1890
1891       if(tnh) then
1892        HNose1=Hnose_nh(EK,potE)
1893        H0=HNose1
1894        write (iout,*) 'H0= ',H0
1895       endif
1896
1897       if (hmc.gt.0) then
1898          hmc_acc=0
1899          hmc_etot=potE+EK
1900           if(me.eq.king.or..not.out1file)
1901      &       write(iout,*) 'HMC',hmc_etot,potE,EK
1902          do i=1,2*nres
1903            do j=1,3
1904             dc_hmc(j,i)=dc(j,i)
1905            enddo
1906          enddo
1907       endif
1908
1909       call cartgrad
1910       call lagrangian
1911       call max_accel
1912       if (amax*d_time .gt. dvmax) then
1913         d_time=d_time*dvmax/amax
1914         if(me.eq.king.or..not.out1file) write (iout,*) 
1915      &   "Time step reduced to",d_time,
1916      &   " because of too large initial acceleration."
1917       endif
1918       if(me.eq.king.or..not.out1file)then 
1919        write(iout,*) "Potential energy and its components"
1920        call enerprint(potEcomp)
1921 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1922       endif
1923       potE=potEcomp(0)-potEcomp(20)
1924       totE=EK+potE
1925       itime=0
1926       if (ntwe.ne.0) call statout(itime)
1927       if(me.eq.king.or..not.out1file)
1928      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1929      &   " Kinetic energy",EK," potential energy",potE, 
1930      &   " total energy",totE," maximum acceleration ",
1931      &   amax
1932       if (large) then
1933         write (iout,*) "Initial coordinates"
1934         do i=1,nres
1935           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1936      &    (c(j,i+nres),j=1,3)
1937         enddo
1938         write (iout,*) "Initial dC"
1939         do i=0,nres
1940           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1941      &    (dc(j,i+nres),j=1,3)
1942         enddo
1943         write (iout,*) "Initial velocities"
1944         write (iout,"(13x,' backbone ',23x,' side chain')")
1945         do i=0,nres
1946           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1947      &    (d_t(j,i+nres),j=1,3)
1948         enddo
1949         write (iout,*) "Initial accelerations"
1950         do i=0,nres
1951 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1952           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1953      &    (d_a(j,i+nres),j=1,3)
1954         enddo
1955       endif
1956       do i=0,2*nres
1957         do j=1,3
1958           dc_old(j,i)=dc(j,i)
1959           d_t_old(j,i)=d_t(j,i)
1960           d_a_old(j,i)=d_a(j,i)
1961         enddo
1962 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1963       enddo 
1964       if (RESPA) then
1965 #ifdef MPI
1966       tt0 =MPI_Wtime()
1967 #else
1968       tt0 = tcpu()
1969 #endif
1970         call zerograd
1971         call etotal_short(energia_short)
1972 #ifdef TIMING_ENE
1973 #ifdef MPI
1974         t_eshort=t_eshort+MPI_Wtime()-tt0
1975 #else
1976         t_eshort=t_eshort+tcpu()-tt0
1977 #endif
1978 #endif
1979
1980         if(tnp .or. tnp1) then
1981          E_short=energia_short(0)
1982          HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1983          Csplit=Hnose1
1984 c         Csplit =110
1985 c_new_var_csplit          Csplit=H0-E_long 
1986 c          Csplit = H0-energia_short(0)
1987           write(iout,*) 'Csplit= ',Csplit
1988         endif
1989
1990
1991         call cartgrad
1992         call lagrangian
1993         if(.not.out1file .and. large) then
1994           write (iout,*) "energia_long",energia_long(0),
1995      &     " energia_short",energia_short(0),
1996      &     " total",energia_long(0)+energia_short(0)
1997           write (iout,*) "Initial fast-force accelerations"
1998           do i=0,nres
1999             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2000      &      (d_a(j,i+nres),j=1,3)
2001           enddo
2002         endif
2003 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2004         do i=0,2*nres
2005           do j=1,3
2006             d_a_short(j,i)=d_a(j,i)
2007           enddo
2008         enddo
2009 #ifdef MPI
2010         tt0=MPI_Wtime()
2011 #else
2012         tt0=tcpu()
2013 #endif
2014         call zerograd
2015         call etotal_long(energia_long)
2016 #ifdef TIMING_ENE
2017 #ifdef MPI
2018         t_elong=t_elong+MPI_Wtime()-tt0
2019 #else
2020         t_elong=t_elong+tcpu()-tt0
2021 #endif
2022 #endif
2023         call cartgrad
2024         call lagrangian
2025         if(.not.out1file .and. large) then
2026           write (iout,*) "energia_long",energia_long(0)
2027           write (iout,*) "Initial slow-force accelerations"
2028           do i=0,nres
2029             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2030      &      (d_a(j,i+nres),j=1,3)
2031           enddo
2032         endif
2033 #ifdef MPI
2034         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2035 #else
2036         t_enegrad=t_enegrad+tcpu()-tt0
2037 #endif
2038       endif
2039
2040
2041
2042       return
2043       end
2044 c-----------------------------------------------------------
2045       subroutine random_vel
2046       implicit real*8 (a-h,o-z)
2047       include 'DIMENSIONS'
2048       include 'COMMON.CONTROL'
2049       include 'COMMON.VAR'
2050       include 'COMMON.MD'
2051 #ifndef LANG0
2052       include 'COMMON.LANGEVIN'
2053 #else
2054       include 'COMMON.LANGEVIN.lang0'
2055 #endif
2056       include 'COMMON.CHAIN'
2057       include 'COMMON.DERIV'
2058       include 'COMMON.GEO'
2059       include 'COMMON.LOCAL'
2060       include 'COMMON.INTERACT'
2061       include 'COMMON.IOUNITS'
2062       include 'COMMON.NAMES'
2063       include 'COMMON.TIME1'
2064       double precision xv,sigv,lowb,highb
2065 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m 
2066 c First generate velocities in the eigenspace of the G matrix
2067 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2068 c      call flush(iout)
2069 c      write (iout,*) "RANDOM_VEL dimen",dimen
2070       xv=0.0d0
2071       ii=0
2072       do i=1,dimen
2073         do k=1,3
2074           ii=ii+1
2075           sigv=dsqrt((Rb*t_bath)/geigen(i))
2076           lowb=-5*sigv
2077           highb=5*sigv
2078           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2079 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2080 c     &      " d_t_work_new",d_t_work_new(ii)
2081         enddo
2082       enddo
2083       call flush(iout)
2084 c diagnostics
2085 c      Ek1=0.0d0
2086 c      ii=0
2087 c      do i=1,dimen
2088 c        do k=1,3
2089 c          ii=ii+1
2090 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2091 c        enddo
2092 c      enddo
2093 c      write (iout,*) "Ek from eigenvectors",Ek1
2094 c end diagnostics
2095 c Transform velocities to UNRES coordinate space
2096       do k=0,2       
2097         do i=1,dimen
2098           ind=(i-1)*3+k+1
2099           d_t_work(ind)=0.0d0
2100           do j=1,dimen
2101             d_t_work(ind)=d_t_work(ind)
2102      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2103           enddo
2104 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2105 c          call flush(iout)
2106         enddo
2107       enddo
2108 c Transfer to the d_t vector
2109       do j=1,3
2110         d_t(j,0)=d_t_work(j)
2111       enddo 
2112       ind=3
2113       do i=nnt,nct-1
2114         do j=1,3 
2115           ind=ind+1
2116           d_t(j,i)=d_t_work(ind)
2117         enddo
2118       enddo
2119 c      do i=0,nres-1
2120 c        write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2121 c      enddo
2122 c      call flush(iout)
2123       do i=nnt,nct
2124         if (itype(i).ne.10) then
2125           do j=1,3
2126             ind=ind+1
2127             d_t(j,i+nres)=d_t_work(ind)
2128           enddo
2129         endif
2130       enddo
2131 c      call kinetic(EK)
2132 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2133 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2134 c      call flush(iout)
2135       return
2136       end
2137 #ifndef LANG0
2138 c-----------------------------------------------------------
2139       subroutine sd_verlet_p_setup
2140 c Sets up the parameters of stochastic Verlet algorithm       
2141       implicit real*8 (a-h,o-z)
2142       include 'DIMENSIONS'
2143 #ifdef MPI
2144       include 'mpif.h'
2145 #endif
2146       include 'COMMON.CONTROL'
2147       include 'COMMON.VAR'
2148       include 'COMMON.MD'
2149 #ifndef LANG0
2150       include 'COMMON.LANGEVIN'
2151 #else
2152       include 'COMMON.LANGEVIN.lang0'
2153 #endif
2154       include 'COMMON.CHAIN'
2155       include 'COMMON.DERIV'
2156       include 'COMMON.GEO'
2157       include 'COMMON.LOCAL'
2158       include 'COMMON.INTERACT'
2159       include 'COMMON.IOUNITS'
2160       include 'COMMON.NAMES'
2161       include 'COMMON.TIME1'
2162       double precision emgdt(MAXRES6),
2163      & pterm,vterm,rho,rhoc,vsig,
2164      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2165      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2166      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2167       logical lprn /.false./
2168       double precision zero /1.0d-8/, gdt_radius /0.05d0/ 
2169       double precision ktm
2170 #ifdef MPI
2171       tt0 = MPI_Wtime()
2172 #else
2173       tt0 = tcpu()
2174 #endif
2175 c
2176 c AL 8/17/04 Code adapted from tinker
2177 c
2178 c Get the frictional and random terms for stochastic dynamics in the
2179 c eigenspace of mass-scaled UNRES friction matrix
2180 c
2181       do i = 1, dimen
2182             gdt = fricgam(i) * d_time
2183 c
2184 c Stochastic dynamics reduces to simple MD for zero friction
2185 c
2186             if (gdt .le. zero) then
2187                pfric_vec(i) = 1.0d0
2188                vfric_vec(i) = d_time
2189                afric_vec(i) = 0.5d0 * d_time * d_time
2190                prand_vec(i) = 0.0d0
2191                vrand_vec1(i) = 0.0d0
2192                vrand_vec2(i) = 0.0d0
2193 c
2194 c Analytical expressions when friction coefficient is large
2195 c
2196             else 
2197                if (gdt .ge. gdt_radius) then
2198                   egdt = dexp(-gdt)
2199                   pfric_vec(i) = egdt
2200                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2201                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2202                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2203                   vterm = 1.0d0 - egdt**2
2204                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2205 c
2206 c Use series expansions when friction coefficient is small
2207 c
2208                else
2209                   gdt2 = gdt * gdt
2210                   gdt3 = gdt * gdt2
2211                   gdt4 = gdt2 * gdt2
2212                   gdt5 = gdt2 * gdt3
2213                   gdt6 = gdt3 * gdt3
2214                   gdt7 = gdt3 * gdt4
2215                   gdt8 = gdt4 * gdt4
2216                   gdt9 = gdt4 * gdt5
2217                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2218      &                          - gdt5/120.0d0 + gdt6/720.0d0
2219      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2220      &                          - gdt9/362880.0d0) / fricgam(i)**2
2221                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2222                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2223                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2224      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2225      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2226      &                       + 127.0d0*gdt9/90720.0d0
2227                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2228      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2229      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2230      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2231                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2232      &                       - 17.0d0*gdt2/1280.0d0
2233      &                       + 17.0d0*gdt3/6144.0d0
2234      &                       + 40967.0d0*gdt4/34406400.0d0
2235      &                       - 57203.0d0*gdt5/275251200.0d0
2236      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2237                end if
2238 c
2239 c Compute the scaling factors of random terms for the nonzero friction case
2240 c
2241                ktm = 0.5d0*d_time/fricgam(i)
2242                psig = dsqrt(ktm*pterm) / fricgam(i)
2243                vsig = dsqrt(ktm*vterm)
2244                rhoc = dsqrt(1.0d0 - rho*rho)
2245                prand_vec(i) = psig 
2246                vrand_vec1(i) = vsig * rho 
2247                vrand_vec2(i) = vsig * rhoc
2248             end if
2249       end do
2250       if (lprn) then
2251       write (iout,*) 
2252      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2253      &  " vrand_vec2"
2254       do i=1,dimen
2255         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2256      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2257       enddo
2258       endif
2259 c
2260 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2261 c
2262 #ifndef   LANG0
2263       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2264       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2265       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2266       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2267       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2268       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2269 #endif
2270 #ifdef MPI
2271       t_sdsetup=t_sdsetup+MPI_Wtime()
2272 #else
2273       t_sdsetup=t_sdsetup+tcpu()-tt0
2274 #endif
2275       return
2276       end
2277 c-------------------------------------------------------------      
2278       subroutine eigtransf1(n,ndim,ab,d,c)
2279       implicit none
2280       integer n,ndim
2281       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2282       integer i,j,k
2283       do i=1,n
2284         do j=1,n
2285           c(i,j)=0.0d0
2286           do k=1,n
2287             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2288           enddo
2289         enddo
2290       enddo
2291       return
2292       end
2293 c-------------------------------------------------------------      
2294       subroutine eigtransf(n,ndim,a,b,d,c)
2295       implicit none
2296       integer n,ndim
2297       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2298       integer i,j,k
2299       do i=1,n
2300         do j=1,n
2301           c(i,j)=0.0d0
2302           do k=1,n
2303             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2304           enddo
2305         enddo
2306       enddo
2307       return
2308       end
2309 c-------------------------------------------------------------      
2310       subroutine sd_verlet1
2311 c Applying stochastic velocity Verlet algorithm - step 1 to velocities        
2312       implicit real*8 (a-h,o-z)
2313       include 'DIMENSIONS'
2314       include 'COMMON.CONTROL'
2315       include 'COMMON.VAR'
2316       include 'COMMON.MD'
2317 #ifndef LANG0
2318       include 'COMMON.LANGEVIN'
2319 #else
2320       include 'COMMON.LANGEVIN.lang0'
2321 #endif
2322       include 'COMMON.CHAIN'
2323       include 'COMMON.DERIV'
2324       include 'COMMON.GEO'
2325       include 'COMMON.LOCAL'
2326       include 'COMMON.INTERACT'
2327       include 'COMMON.IOUNITS'
2328       include 'COMMON.NAMES'
2329       double precision stochforcvec(MAXRES6)
2330       common /stochcalc/ stochforcvec
2331       logical lprn /.false./
2332
2333 c      write (iout,*) "dc_old"
2334 c      do i=0,nres
2335 c        write (iout,'(i5,3f10.5,5x,3f10.5)') 
2336 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2337 c      enddo
2338       do j=1,3
2339         dc_work(j)=dc_old(j,0)
2340         d_t_work(j)=d_t_old(j,0)
2341         d_a_work(j)=d_a_old(j,0)
2342       enddo
2343       ind=3
2344       do i=nnt,nct-1
2345         do j=1,3
2346           dc_work(ind+j)=dc_old(j,i)
2347           d_t_work(ind+j)=d_t_old(j,i)
2348           d_a_work(ind+j)=d_a_old(j,i)
2349         enddo
2350         ind=ind+3
2351       enddo
2352       do i=nnt,nct
2353         if (itype(i).ne.10) then
2354           do j=1,3
2355             dc_work(ind+j)=dc_old(j,i+nres)
2356             d_t_work(ind+j)=d_t_old(j,i+nres)
2357             d_a_work(ind+j)=d_a_old(j,i+nres)
2358           enddo
2359           ind=ind+3
2360         endif
2361       enddo
2362 #ifndef LANG0
2363       if (lprn) then
2364       write (iout,*) 
2365      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2366      &  " vrand_mat2"
2367       do i=1,dimen
2368         do j=1,dimen
2369           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2370      &      vfric_mat(i,j),afric_mat(i,j),
2371      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2372         enddo
2373       enddo
2374       endif
2375       do i=1,dimen
2376         ddt1=0.0d0
2377         ddt2=0.0d0
2378         do j=1,dimen
2379           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2380      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2381           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2382           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2383         enddo
2384         d_t_work_new(i)=ddt1+0.5d0*ddt2
2385         d_t_work(i)=ddt1+ddt2
2386       enddo
2387 #endif
2388       do j=1,3
2389         dc(j,0)=dc_work(j)
2390         d_t(j,0)=d_t_work(j)
2391       enddo
2392       ind=3     
2393       do i=nnt,nct-1    
2394         do j=1,3
2395           dc(j,i)=dc_work(ind+j)
2396           d_t(j,i)=d_t_work(ind+j)
2397         enddo
2398         ind=ind+3
2399       enddo
2400       do i=nnt,nct
2401         if (itype(i).ne.10) then
2402           inres=i+nres
2403           do j=1,3
2404             dc(j,inres)=dc_work(ind+j)
2405             d_t(j,inres)=d_t_work(ind+j)
2406           enddo
2407           ind=ind+3
2408         endif      
2409       enddo 
2410       return
2411       end
2412 c--------------------------------------------------------------------------
2413       subroutine sd_verlet2
2414 c  Calculating the adjusted velocities for accelerations
2415       implicit real*8 (a-h,o-z)
2416       include 'DIMENSIONS'
2417       include 'COMMON.CONTROL'
2418       include 'COMMON.VAR'
2419       include 'COMMON.MD'
2420 #ifndef LANG0
2421       include 'COMMON.LANGEVIN'
2422 #else
2423       include 'COMMON.LANGEVIN.lang0'
2424 #endif
2425       include 'COMMON.CHAIN'
2426       include 'COMMON.DERIV'
2427       include 'COMMON.GEO'
2428       include 'COMMON.LOCAL'
2429       include 'COMMON.INTERACT'
2430       include 'COMMON.IOUNITS'
2431       include 'COMMON.NAMES'
2432       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2433       common /stochcalc/ stochforcvec
2434 c
2435 c Compute the stochastic forces which contribute to velocity change
2436 c
2437       call stochastic_force(stochforcvecV)
2438
2439 #ifndef LANG0
2440       do i=1,dimen
2441         ddt1=0.0d0
2442         ddt2=0.0d0
2443         do j=1,dimen
2444           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2445           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2446      &     vrand_mat2(i,j)*stochforcvecV(j)
2447         enddo
2448         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2449       enddo
2450 #endif
2451       do j=1,3
2452         d_t(j,0)=d_t_work(j)
2453       enddo
2454       ind=3
2455       do i=nnt,nct-1
2456         do j=1,3
2457           d_t(j,i)=d_t_work(ind+j)
2458         enddo
2459         ind=ind+3
2460       enddo
2461       do i=nnt,nct
2462         if (itype(i).ne.10) then
2463           inres=i+nres
2464           do j=1,3
2465             d_t(j,inres)=d_t_work(ind+j)
2466           enddo
2467           ind=ind+3
2468         endif
2469       enddo 
2470       return
2471       end
2472 c-----------------------------------------------------------
2473       subroutine sd_verlet_ciccotti_setup
2474 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's 
2475 c version 
2476       implicit real*8 (a-h,o-z)
2477       include 'DIMENSIONS'
2478 #ifdef MPI
2479       include 'mpif.h'
2480 #endif
2481       include 'COMMON.CONTROL'
2482       include 'COMMON.VAR'
2483       include 'COMMON.MD'
2484 #ifndef LANG0
2485       include 'COMMON.LANGEVIN'
2486 #else
2487       include 'COMMON.LANGEVIN.lang0'
2488 #endif
2489       include 'COMMON.CHAIN'
2490       include 'COMMON.DERIV'
2491       include 'COMMON.GEO'
2492       include 'COMMON.LOCAL'
2493       include 'COMMON.INTERACT'
2494       include 'COMMON.IOUNITS'
2495       include 'COMMON.NAMES'
2496       include 'COMMON.TIME1'
2497       double precision emgdt(MAXRES6),
2498      & pterm,vterm,rho,rhoc,vsig,
2499      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2500      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2501      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2502       logical lprn /.false./
2503       double precision zero /1.0d-8/, gdt_radius /0.05d0/ 
2504       double precision ktm
2505 #ifdef MPI
2506       tt0 = MPI_Wtime()
2507 #else
2508       tt0 = tcpu()
2509 #endif
2510 c
2511 c AL 8/17/04 Code adapted from tinker
2512 c
2513 c Get the frictional and random terms for stochastic dynamics in the
2514 c eigenspace of mass-scaled UNRES friction matrix
2515 c
2516       do i = 1, dimen
2517             write (iout,*) "i",i," fricgam",fricgam(i)
2518             gdt = fricgam(i) * d_time
2519 c
2520 c Stochastic dynamics reduces to simple MD for zero friction
2521 c
2522             if (gdt .le. zero) then
2523                pfric_vec(i) = 1.0d0
2524                vfric_vec(i) = d_time
2525                afric_vec(i) = 0.5d0*d_time*d_time
2526                prand_vec(i) = afric_vec(i)
2527                vrand_vec2(i) = vfric_vec(i)
2528 c
2529 c Analytical expressions when friction coefficient is large
2530 c
2531             else 
2532                egdt = dexp(-gdt)
2533                pfric_vec(i) = egdt
2534                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2535                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2536                prand_vec(i) = afric_vec(i)
2537                vrand_vec2(i) = vfric_vec(i)
2538 c
2539 c Compute the scaling factors of random terms for the nonzero friction case
2540 c
2541 c               ktm = 0.5d0*d_time/fricgam(i)
2542 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2543 c               vsig = dsqrt(ktm*vterm)
2544 c               prand_vec(i) = psig*afric_vec(i) 
2545 c               vrand_vec2(i) = vsig*vfric_vec(i)
2546             end if
2547       end do
2548       if (lprn) then
2549       write (iout,*) 
2550      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2551      &  " vrand_vec2"
2552       do i=1,dimen
2553         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2554      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2555       enddo
2556       endif
2557 c
2558 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2559 c
2560       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2561       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2562       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2563       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2564       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2565 #ifdef MPI
2566       t_sdsetup=t_sdsetup+MPI_Wtime()
2567 #else
2568       t_sdsetup=t_sdsetup+tcpu()-tt0
2569 #endif
2570       return
2571       end
2572 c-------------------------------------------------------------      
2573       subroutine sd_verlet1_ciccotti
2574 c Applying stochastic velocity Verlet algorithm - step 1 to velocities        
2575       implicit real*8 (a-h,o-z)
2576       include 'DIMENSIONS'
2577 #ifdef MPI
2578       include 'mpif.h'
2579 #endif
2580       include 'COMMON.CONTROL'
2581       include 'COMMON.VAR'
2582       include 'COMMON.MD'
2583 #ifndef LANG0
2584       include 'COMMON.LANGEVIN'
2585 #else
2586       include 'COMMON.LANGEVIN.lang0'
2587 #endif
2588       include 'COMMON.CHAIN'
2589       include 'COMMON.DERIV'
2590       include 'COMMON.GEO'
2591       include 'COMMON.LOCAL'
2592       include 'COMMON.INTERACT'
2593       include 'COMMON.IOUNITS'
2594       include 'COMMON.NAMES'
2595       double precision stochforcvec(MAXRES6)
2596       common /stochcalc/ stochforcvec
2597       logical lprn /.false./
2598
2599 c      write (iout,*) "dc_old"
2600 c      do i=0,nres
2601 c        write (iout,'(i5,3f10.5,5x,3f10.5)') 
2602 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2603 c      enddo
2604       do j=1,3
2605         dc_work(j)=dc_old(j,0)
2606         d_t_work(j)=d_t_old(j,0)
2607         d_a_work(j)=d_a_old(j,0)
2608       enddo
2609       ind=3
2610       do i=nnt,nct-1
2611         do j=1,3
2612           dc_work(ind+j)=dc_old(j,i)
2613           d_t_work(ind+j)=d_t_old(j,i)
2614           d_a_work(ind+j)=d_a_old(j,i)
2615         enddo
2616         ind=ind+3
2617       enddo
2618       do i=nnt,nct
2619         if (itype(i).ne.10) then
2620           do j=1,3
2621             dc_work(ind+j)=dc_old(j,i+nres)
2622             d_t_work(ind+j)=d_t_old(j,i+nres)
2623             d_a_work(ind+j)=d_a_old(j,i+nres)
2624           enddo
2625           ind=ind+3
2626         endif
2627       enddo
2628
2629 #ifndef LANG0
2630       if (lprn) then
2631       write (iout,*) 
2632      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2633      &  " vrand_mat2"
2634       do i=1,dimen
2635         do j=1,dimen
2636           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2637      &      vfric_mat(i,j),afric_mat(i,j),
2638      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2639         enddo
2640       enddo
2641       endif
2642       do i=1,dimen
2643         ddt1=0.0d0
2644         ddt2=0.0d0
2645         do j=1,dimen
2646           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2647      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2648           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2649           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2650         enddo
2651         d_t_work_new(i)=ddt1+0.5d0*ddt2
2652         d_t_work(i)=ddt1+ddt2
2653       enddo
2654 #endif
2655       do j=1,3
2656         dc(j,0)=dc_work(j)
2657         d_t(j,0)=d_t_work(j)
2658       enddo
2659       ind=3     
2660       do i=nnt,nct-1    
2661         do j=1,3
2662           dc(j,i)=dc_work(ind+j)
2663           d_t(j,i)=d_t_work(ind+j)
2664         enddo
2665         ind=ind+3
2666       enddo
2667       do i=nnt,nct
2668         if (itype(i).ne.10) then
2669           inres=i+nres
2670           do j=1,3
2671             dc(j,inres)=dc_work(ind+j)
2672             d_t(j,inres)=d_t_work(ind+j)
2673           enddo
2674           ind=ind+3
2675         endif      
2676       enddo 
2677       return
2678       end
2679 c--------------------------------------------------------------------------
2680       subroutine sd_verlet2_ciccotti
2681 c  Calculating the adjusted velocities for accelerations
2682       implicit real*8 (a-h,o-z)
2683       include 'DIMENSIONS'
2684       include 'COMMON.CONTROL'
2685       include 'COMMON.VAR'
2686       include 'COMMON.MD'
2687 #ifndef LANG0
2688       include 'COMMON.LANGEVIN'
2689 #else
2690       include 'COMMON.LANGEVIN.lang0'
2691 #endif
2692       include 'COMMON.CHAIN'
2693       include 'COMMON.DERIV'
2694       include 'COMMON.GEO'
2695       include 'COMMON.LOCAL'
2696       include 'COMMON.INTERACT'
2697       include 'COMMON.IOUNITS'
2698       include 'COMMON.NAMES'
2699       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2700       common /stochcalc/ stochforcvec
2701 c
2702 c Compute the stochastic forces which contribute to velocity change
2703 c
2704       call stochastic_force(stochforcvecV)
2705 #ifndef LANG0
2706       do i=1,dimen
2707         ddt1=0.0d0
2708         ddt2=0.0d0
2709         do j=1,dimen
2710
2711           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2712 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2713           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2714         enddo
2715         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2716       enddo
2717 #endif
2718       do j=1,3
2719         d_t(j,0)=d_t_work(j)
2720       enddo
2721       ind=3
2722       do i=nnt,nct-1
2723         do j=1,3
2724           d_t(j,i)=d_t_work(ind+j)
2725         enddo
2726         ind=ind+3
2727       enddo
2728       do i=nnt,nct
2729         if (itype(i).ne.10) then
2730           inres=i+nres
2731           do j=1,3
2732             d_t(j,inres)=d_t_work(ind+j)
2733           enddo
2734           ind=ind+3
2735         endif
2736       enddo 
2737       return
2738       end
2739 #endif
2740 c------------------------------------------------------
2741       double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2742       implicit none
2743       double precision ek,s,e,pi,Q,t_bath,Rb
2744       integer dimenl
2745       Rb=0.001986d0
2746       HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2747 c      print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2748 c     &      pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2749       return
2750       end
2751 c-----------------------------------------------------------------
2752       double precision function HNose_nh(eki,e)
2753       implicit real*8 (a-h,o-z)
2754       include 'DIMENSIONS'
2755       include 'COMMON.MD'
2756       HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2757       do i=2,nnos
2758         HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2759       enddo
2760 c      write(4,'(5e15.5)') 
2761 c     &       vlogs(1),xlogs(1),HNose,eki,e
2762       return
2763       end
2764 c-----------------------------------------------------------------
2765       SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2766       implicit real*8 (a-h,o-z)
2767       include 'DIMENSIONS'
2768       include 'COMMON.MD'
2769       double precision akin,gnkt,dt,aa,gkt,scale
2770       double precision wdti(maxyosh),wdti2(maxyosh),
2771      &                 wdti4(maxyosh),wdti8(maxyosh)
2772       integer i,iresn,iyosh,inos,nnos1
2773
2774       dt=d_time
2775       nnos1=nnos+1
2776       GKT = Rb*t_bath
2777       GNKT = dimen3*GKT
2778       akin=akin*2
2779
2780       
2781 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2782 C INTEGRATION FROM t=0 TO t=DT/2
2783 C GET THE TOTAL KINETIC ENERGY
2784       SCALE = 1.D0
2785 c      CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2786 C UPDATE THE FORCES
2787       GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2788 C START THE MULTIPLE TIME STEP PROCEDURE
2789       DO IRESN = 1,NRESN
2790        DO IYOSH = 1,NYOSH
2791 C UPDATE THE THERMOSTAT VELOCITIES
2792         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2793         DO INOS = 1,NNOS-1
2794          AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2795          VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2796      &          + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2797         ENDDO
2798 C UPDATE THE PARTICLE VELOCITIES
2799         AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2800         SCALE = SCALE*AA
2801 C UPDATE THE FORCES
2802         GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2803 C UPDATE THE THERMOSTAT POSITIONS
2804         DO INOS = 1,NNOS
2805          XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2806         ENDDO
2807 C UPDATE THE THERMOSTAT VELOCITIES
2808         DO INOS = 1,NNOS-1
2809          AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2810          VLOGS(INOS) = VLOGS(INOS)*AA*AA
2811      &      + WDTI4(IYOSH)*GLOGS(INOS)*AA
2812          GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2813      &      -GKT)/QMASS(INOS+1)
2814         ENDDO
2815         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2816        ENDDO
2817       ENDDO
2818 C UPDATE THE PARTICLE VELOCITIES
2819 c outside of this subroutine
2820 c      DO I = 1,N
2821 c       VX(I) = VX(I)*SCALE
2822 c       VY(I) = VY(I)*SCALE
2823 c       VZ(I) = VZ(I)*SCALE
2824 c      ENDDO
2825       RETURN
2826       END
2827 c-----------------------------------------------------------------
2828       subroutine tnp1_respa_i_step1
2829 c Applying Nose-Poincare algorithm - step 1 to coordinates
2830 c JPSJ 70 75 (2001) S. Nose
2831 c
2832 c d_t is not updated here
2833 c
2834       implicit real*8 (a-h,o-z)
2835       include 'DIMENSIONS'
2836       include 'COMMON.CONTROL'
2837       include 'COMMON.VAR'
2838       include 'COMMON.MD'
2839       include 'COMMON.CHAIN'
2840       include 'COMMON.DERIV'
2841       include 'COMMON.GEO'
2842       include 'COMMON.LOCAL'
2843       include 'COMMON.INTERACT'
2844       include 'COMMON.IOUNITS'
2845       include 'COMMON.NAMES'
2846       double precision adt,adt2,tmp
2847         
2848       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2849       s12_np=s_np*tmp**2
2850       pistar=pi_np/tmp
2851       s12_dt=d_time/s12_np
2852       d_time_s12=d_time*0.5*s12_np
2853
2854       do j=1,3
2855         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2856         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2857       enddo
2858       do i=nnt,nct-1    
2859         do j=1,3    
2860           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2861           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2862         enddo
2863       enddo
2864       do i=nnt,nct
2865         if (itype(i).ne.10) then
2866           inres=i+nres
2867           do j=1,3    
2868            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2869            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2870           enddo
2871         endif      
2872       enddo 
2873       return
2874       end
2875 c---------------------------------------------------------------------
2876       subroutine tnp1_respa_i_step2
2877 c  Step 2 of the velocity Verlet algorithm: update velocities
2878       implicit real*8 (a-h,o-z)
2879       include 'DIMENSIONS'
2880       include 'COMMON.CONTROL'
2881       include 'COMMON.VAR'
2882       include 'COMMON.MD'
2883       include 'COMMON.CHAIN'
2884       include 'COMMON.DERIV'
2885       include 'COMMON.GEO'
2886       include 'COMMON.LOCAL'
2887       include 'COMMON.INTERACT'
2888       include 'COMMON.IOUNITS'
2889       include 'COMMON.NAMES'
2890
2891       double precision d_time_s12
2892
2893       do i=0,2*nres
2894        do j=1,3
2895         d_t(j,i)=d_t_new(j,i)
2896        enddo
2897       enddo
2898
2899       call kinetic(EK)
2900       EK=EK/s12_np**2
2901
2902       d_time_s12=0.5d0*s12_np*d_time
2903
2904       do j=1,3
2905         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2906       enddo
2907       do i=nnt,nct-1
2908         do j=1,3
2909           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2910         enddo
2911       enddo
2912       do i=nnt,nct
2913         if (itype(i).ne.10) then
2914           inres=i+nres
2915           do j=1,3
2916             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2917           enddo
2918         endif
2919       enddo 
2920
2921       pistar=pistar+(EK-0.5*(E_old+potE)
2922      &       -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2923       tmp=1+pistar/(2*Q_np)*0.5*d_time
2924       s_np=s12_np*tmp**2
2925       pi_np=pistar/tmp
2926
2927       return
2928       end
2929 c-------------------------------------------------------
2930
2931       subroutine tnp1_step1
2932 c Applying Nose-Poincare algorithm - step 1 to coordinates
2933 c JPSJ 70 75 (2001) S. Nose
2934 c
2935 c d_t is not updated here
2936 c
2937       implicit real*8 (a-h,o-z)
2938       include 'DIMENSIONS'
2939       include 'COMMON.CONTROL'
2940       include 'COMMON.VAR'
2941       include 'COMMON.MD'
2942       include 'COMMON.CHAIN'
2943       include 'COMMON.DERIV'
2944       include 'COMMON.GEO'
2945       include 'COMMON.LOCAL'
2946       include 'COMMON.INTERACT'
2947       include 'COMMON.IOUNITS'
2948       include 'COMMON.NAMES'
2949       double precision adt,adt2,tmp
2950         
2951       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2952       s12_np=s_np*tmp**2
2953       pistar=pi_np/tmp
2954       s12_dt=d_time/s12_np
2955       d_time_s12=d_time*0.5*s12_np
2956
2957       do j=1,3
2958         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2959         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2960       enddo
2961       do i=nnt,nct-1    
2962         do j=1,3    
2963           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2964           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2965         enddo
2966       enddo
2967       do i=nnt,nct
2968         if (itype(i).ne.10) then
2969           inres=i+nres
2970           do j=1,3    
2971            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2972            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2973           enddo
2974         endif      
2975       enddo 
2976       return
2977       end
2978 c---------------------------------------------------------------------
2979       subroutine tnp1_step2
2980 c  Step 2 of the velocity Verlet algorithm: update velocities
2981       implicit real*8 (a-h,o-z)
2982       include 'DIMENSIONS'
2983       include 'COMMON.CONTROL'
2984       include 'COMMON.VAR'
2985       include 'COMMON.MD'
2986       include 'COMMON.CHAIN'
2987       include 'COMMON.DERIV'
2988       include 'COMMON.GEO'
2989       include 'COMMON.LOCAL'
2990       include 'COMMON.INTERACT'
2991       include 'COMMON.IOUNITS'
2992       include 'COMMON.NAMES'
2993
2994       double precision d_time_s12
2995
2996       do i=0,2*nres
2997        do j=1,3
2998         d_t(j,i)=d_t_new(j,i)
2999        enddo
3000       enddo
3001
3002       call kinetic(EK)
3003       EK=EK/s12_np**2
3004
3005       d_time_s12=0.5d0*s12_np*d_time
3006
3007       do j=1,3
3008         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3009       enddo
3010       do i=nnt,nct-1
3011         do j=1,3
3012           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3013         enddo
3014       enddo
3015       do i=nnt,nct
3016         if (itype(i).ne.10) then
3017           inres=i+nres
3018           do j=1,3
3019             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3020           enddo
3021         endif
3022       enddo 
3023
3024 cd      write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3025       pistar=pistar+(EK-0.5*(E_old+potE)
3026      &       -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3027       tmp=1+pistar/(2*Q_np)*0.5*d_time
3028       s_np=s12_np*tmp**2
3029       pi_np=pistar/tmp
3030
3031       return
3032       end
3033
3034 c-----------------------------------------------------------------
3035       subroutine tnp_respa_i_step1
3036 c Applying Nose-Poincare algorithm - step 1 to coordinates
3037 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3038 c
3039 c d_t is not updated here, it is destroyed
3040 c
3041       implicit real*8 (a-h,o-z)
3042       include 'DIMENSIONS'
3043       include 'COMMON.CONTROL'
3044       include 'COMMON.VAR'
3045       include 'COMMON.MD'
3046       include 'COMMON.CHAIN'
3047       include 'COMMON.DERIV'
3048       include 'COMMON.GEO'
3049       include 'COMMON.LOCAL'
3050       include 'COMMON.INTERACT'
3051       include 'COMMON.IOUNITS'
3052       include 'COMMON.NAMES'
3053       double precision C_np,d_time_s,tmp,d_time_ss
3054
3055       d_time_s=d_time*0.5*s_np        
3056 ct2      d_time_s=d_time*0.5*s12_np
3057
3058       do j=1,3
3059         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3060       enddo
3061       do i=nnt,nct-1    
3062         do j=1,3    
3063           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3064         enddo
3065       enddo
3066       do i=nnt,nct
3067         if (itype(i).ne.10) then
3068           inres=i+nres
3069           do j=1,3    
3070            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3071           enddo
3072         endif      
3073       enddo 
3074
3075       do i=0,2*nres
3076        do j=1,3
3077         d_t(j,i)=d_t_new(j,i)
3078        enddo
3079       enddo
3080
3081       call kinetic(EK)
3082       EK=EK/s_np**2
3083
3084       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3085      &                     -pi_np
3086
3087       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3088       tmp=0.5*d_time*pistar/Q_np
3089       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3090
3091       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3092 ct2      d_time_ss=d_time/s12_np
3093 c      d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np) 
3094
3095       do j=1,3
3096         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3097       enddo
3098       do i=nnt,nct-1    
3099         do j=1,3    
3100           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3101         enddo
3102       enddo
3103       do i=nnt,nct
3104         if (itype(i).ne.10) then
3105           inres=i+nres
3106           do j=1,3    
3107            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3108           enddo
3109         endif      
3110       enddo 
3111
3112       return
3113       end
3114 c---------------------------------------------------------------------
3115
3116       subroutine tnp_respa_i_step2
3117 c  Step 2 of the velocity Verlet algorithm: update velocities
3118       implicit real*8 (a-h,o-z)
3119       include 'DIMENSIONS'
3120       include 'COMMON.CONTROL'
3121       include 'COMMON.VAR'
3122       include 'COMMON.MD'
3123       include 'COMMON.CHAIN'
3124       include 'COMMON.DERIV'
3125       include 'COMMON.GEO'
3126       include 'COMMON.LOCAL'
3127       include 'COMMON.INTERACT'
3128       include 'COMMON.IOUNITS'
3129       include 'COMMON.NAMES'
3130
3131       double precision d_time_s
3132
3133       EK=EK*(s_np/s12_np)**2
3134       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3135       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3136      &                              -HNose1+Csplit)         
3137
3138 cr      print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3139       d_time_s=d_time*0.5*s12_np
3140 c      d_time_s=d_time*0.5*s_np
3141
3142       do j=1,3
3143         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3144       enddo
3145       do i=nnt,nct-1
3146         do j=1,3
3147           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3148         enddo
3149       enddo
3150       do i=nnt,nct
3151         if (itype(i).ne.10) then
3152           inres=i+nres
3153           do j=1,3
3154             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3155           enddo
3156         endif
3157       enddo 
3158
3159       s_np=s12_np
3160
3161       return
3162       end
3163 c-----------------------------------------------------------------
3164       subroutine tnp_respa_step1
3165 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3166 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3167 c
3168 c d_t is not updated here, it is destroyed
3169 c
3170       implicit real*8 (a-h,o-z)
3171       include 'DIMENSIONS'
3172       include 'COMMON.CONTROL'
3173       include 'COMMON.VAR'
3174       include 'COMMON.MD'
3175       include 'COMMON.CHAIN'
3176       include 'COMMON.DERIV'
3177       include 'COMMON.GEO'
3178       include 'COMMON.LOCAL'
3179       include 'COMMON.INTERACT'
3180       include 'COMMON.IOUNITS'
3181       include 'COMMON.NAMES'
3182       double precision C_np,d_time_s,tmp,d_time_ss
3183       double precision energia(0:n_ene)
3184
3185       d_time_s=d_time*0.5*s_np        
3186
3187       do j=1,3
3188         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3189       enddo
3190       do i=nnt,nct-1    
3191         do j=1,3    
3192           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3193         enddo
3194       enddo
3195       do i=nnt,nct
3196         if (itype(i).ne.10) then
3197           inres=i+nres
3198           do j=1,3    
3199            d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3200           enddo
3201         endif      
3202       enddo 
3203
3204
3205 c      C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3206 c     &                     -pi_np
3207 c
3208 c      pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3209 c      tmp=0.5*d_time*pistar/Q_np
3210 c      s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3211 c      write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3212
3213 ct1      pi_np=pistar
3214 c      sold_np=s_np
3215 c      s_np=s12_np
3216
3217 c-------------------------------------
3218 c test of reviewer's comment
3219        pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3220 cr       print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3221 c-------------------------------------
3222
3223       return
3224       end
3225 c---------------------------------------------------------------------
3226       subroutine tnp_respa_step2
3227 c  Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3228       implicit real*8 (a-h,o-z)
3229       include 'DIMENSIONS'
3230       include 'COMMON.CONTROL'
3231       include 'COMMON.VAR'
3232       include 'COMMON.MD'
3233       include 'COMMON.CHAIN'
3234       include 'COMMON.DERIV'
3235       include 'COMMON.GEO'
3236       include 'COMMON.LOCAL'
3237       include 'COMMON.INTERACT'
3238       include 'COMMON.IOUNITS'
3239       include 'COMMON.NAMES'
3240
3241       double precision d_time_s
3242
3243 ct1      s12_np=s_np
3244 ct2      pistar=pi_np
3245
3246 ct      call kinetic(EK)
3247 ct      HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3248 ct      pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3249 ct     &                              -0.5*d_time*(HNose1-H0)         
3250
3251 c-------------------------------------
3252 c test of reviewer's comment
3253       pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3254 cr      print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3255 c-------------------------------------
3256       d_time_s=d_time*0.5*s_np
3257
3258       do j=1,3
3259         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3260       enddo
3261       do i=nnt,nct-1
3262         do j=1,3
3263           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3264         enddo
3265       enddo
3266       do i=nnt,nct
3267         if (itype(i).ne.10) then
3268           inres=i+nres
3269           do j=1,3
3270             d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3271           enddo
3272         endif
3273       enddo 
3274
3275 cd      s_np=s12_np
3276
3277       return
3278       end
3279 c---------------------------------------------------------------------
3280       subroutine tnp_step1
3281 c Applying Nose-Poincare algorithm - step 1 to coordinates
3282 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3283 c
3284 c d_t is not updated here, it is destroyed
3285 c
3286       implicit real*8 (a-h,o-z)
3287       include 'DIMENSIONS'
3288       include 'COMMON.CONTROL'
3289       include 'COMMON.VAR'
3290       include 'COMMON.MD'
3291       include 'COMMON.CHAIN'
3292       include 'COMMON.DERIV'
3293       include 'COMMON.GEO'
3294       include 'COMMON.LOCAL'
3295       include 'COMMON.INTERACT'
3296       include 'COMMON.IOUNITS'
3297       include 'COMMON.NAMES'
3298       double precision C_np,d_time_s,tmp,d_time_ss
3299
3300       d_time_s=d_time*0.5*s_np        
3301
3302       do j=1,3
3303         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3304       enddo
3305       do i=nnt,nct-1    
3306         do j=1,3    
3307           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3308         enddo
3309       enddo
3310       do i=nnt,nct
3311         if (itype(i).ne.10) then
3312           inres=i+nres
3313           do j=1,3    
3314            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3315           enddo
3316         endif      
3317       enddo 
3318
3319       do i=0,2*nres
3320        do j=1,3
3321         d_t(j,i)=d_t_new(j,i)
3322        enddo
3323       enddo
3324
3325       call kinetic(EK)
3326       EK=EK/s_np**2
3327
3328       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3329      &                     -pi_np
3330
3331       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3332       tmp=0.5*d_time*pistar/Q_np
3333       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3334 c      write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3335
3336       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3337
3338       do j=1,3
3339         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3340       enddo
3341       do i=nnt,nct-1    
3342         do j=1,3    
3343           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3344         enddo
3345       enddo
3346       do i=nnt,nct
3347         if (itype(i).ne.10) then
3348           inres=i+nres
3349           do j=1,3    
3350            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3351           enddo
3352         endif      
3353       enddo 
3354
3355       return
3356       end
3357 c-----------------------------------------------------------------
3358       subroutine tnp_step2
3359 c  Step 2 of the velocity Verlet algorithm: update velocities
3360       implicit real*8 (a-h,o-z)
3361       include 'DIMENSIONS'
3362       include 'COMMON.CONTROL'
3363       include 'COMMON.VAR'
3364       include 'COMMON.MD'
3365       include 'COMMON.CHAIN'
3366       include 'COMMON.DERIV'
3367       include 'COMMON.GEO'
3368       include 'COMMON.LOCAL'
3369       include 'COMMON.INTERACT'
3370       include 'COMMON.IOUNITS'
3371       include 'COMMON.NAMES'
3372
3373       double precision d_time_s
3374
3375       EK=EK*(s_np/s12_np)**2
3376       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3377       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3378      &                              -0.5*d_time*(HNose1-H0)         
3379
3380 cd      write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3381       d_time_s=d_time*0.5*s12_np
3382
3383       do j=1,3
3384         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3385       enddo
3386       do i=nnt,nct-1
3387         do j=1,3
3388           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3389         enddo
3390       enddo
3391       do i=nnt,nct
3392         if (itype(i).ne.10) then
3393           inres=i+nres
3394           do j=1,3
3395             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3396           enddo
3397         endif
3398       enddo 
3399
3400       s_np=s12_np
3401
3402       return
3403       end
3404
3405       subroutine hmc_test(itime)
3406       implicit real*8 (a-h,o-z)
3407       include 'DIMENSIONS'
3408       include 'COMMON.CONTROL'
3409       include 'COMMON.MD'
3410       include 'COMMON.CHAIN'
3411
3412            hmc_acc=hmc_acc+1
3413            delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
3414            if (delta .lt. -50.0d0) then
3415                 delta=0.0d0
3416            else
3417                 delta=dexp(delta)
3418            endif
3419            xxx=ran_number(0.0d0,1.0d0)
3420
3421            if (me.eq.king .or. .not. out1file)
3422      &       write(iout,'(a8,i5,6f10.4)') 
3423      &        'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
3424
3425            if (delta .le. xxx) then
3426             do i=1,2*nres
3427              do j=1,3
3428               dc(j,i)=dc_hmc(j,i)
3429              enddo
3430             enddo
3431             itime=itime-hmc
3432             totT=totThmc
3433            else
3434             if (me.eq.king .or. .not. out1file)
3435      &       write(iout,*) 'HMC accepting new'
3436             totThmc=totT
3437             do i=1,2*nres
3438              do j=1,3
3439               dc_hmc(j,i)=dc(j,i)
3440              enddo
3441             enddo
3442            endif
3443
3444            call chainbuild_cart
3445            call random_vel
3446            do i=0,2*nres
3447             do j=1,3
3448               d_t_old(j,i)=d_t(j,i)
3449             enddo
3450            enddo
3451            call kinetic(EK)
3452            kinetic_T=2.0d0/(dimen3*Rb)*EK
3453            call etotal(potEcomp)
3454            potE=potEcomp(0)
3455            hmc_etot=potE+EK
3456            if (me.eq.king .or. .not. out1file)
3457      &      write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK
3458
3459
3460       return
3461       end