2 c------------------------------------------------
3 c The driver for molecular dynamics subroutines
4 c------------------------------------------------
5 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
10 include 'COMMON.LANGEVIN'
11 include 'COMMON.CHAIN'
12 include 'COMMON.DERIV'
14 include 'COMMON.LOCAL'
15 include 'COMMON.INTERACT'
16 include 'COMMON.IOUNITS'
17 include 'COMMON.NAMES'
18 include 'COMMON.TIME1'
19 double precision cm(3),L(3),vcm(3)
20 double precision energia(0:n_ene)
24 common /gucio/ cm,energia
32 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
34 c Determine the inverse of the inertia matrix.
35 call setup_MD_matrices
38 t_MDsetup = tcpu()-tt0
40 c Entering the MD loop
42 if (lang.eq.2 .or. lang.eq.3) then
45 call sd_verlet_p_setup
47 call sd_verlet_ciccotti_setup
51 pfric0_mat(i,j,0)=pfric_mat(i,j)
52 afric0_mat(i,j,0)=afric_mat(i,j)
53 vfric0_mat(i,j,0)=vfric_mat(i,j)
54 prand0_mat(i,j,0)=prand_mat(i,j)
55 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
56 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
63 else if (lang.eq.1 .or. lang.eq.4) then
66 t_langsetup=tcpu()-tt0
70 if (lang.gt.0 .and. surfarea .and.
71 & mod(itime,reset_fricmat).eq.0) then
72 if (lang.eq.2 .or. lang.eq.3) then
75 call sd_verlet_p_setup
77 call sd_verlet_ciccotti_setup
81 pfric0_mat(i,j,0)=pfric_mat(i,j)
82 afric0_mat(i,j,0)=afric_mat(i,j)
83 vfric0_mat(i,j,0)=vfric_mat(i,j)
84 prand0_mat(i,j,0)=prand_mat(i,j)
85 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
86 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
93 else if (lang.eq.1 .or. lang.eq.4) then
96 write (iout,'(a,i10)')
97 & "Friction matrix reset based on surface area, itime",itime
99 if (reset_vel .and. tbf .and. lang.eq.0
100 & .and. mod(itime,count_reset_vel).eq.0) then
102 write(iout,'(a,f20.2)')
103 & "Velocities reset to random values, time",totT
106 d_t_old(j,i)=d_t(j,i)
110 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
114 d_t(j,0)=d_t(j,0)-vcm(j)
117 kinetic_T=2.0d0/(dimen*Rb)*EK
118 scalfac=dsqrt(T_bath/kinetic_T)
119 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
122 d_t_old(j,i)=scalfac*d_t(j,i)
128 c Time-reversible RESPA algorithm
129 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
130 call RESPA_step(itime)
132 c Variable time step algorithm.
133 call velverlet_step(itime)
136 call brown_step(itime)
138 if (mod(itime,ntwe).eq.0) call statout(itime)
139 if (mod(itime,ntwx).eq.0) then
140 write (tytul,'("time",f8.2)') totT
142 call pdbout(potE,tytul,ipdb)
147 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
148 open(irest2,file=rest2name,status='unknown')
149 write(irest2,*) totT,EK,potE,totE
151 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
154 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
161 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
163 & 'MD calculations setup:',t_MDsetup,
164 & 'Energy & gradient evaluation:',t_enegrad,
165 & 'Stochastic MD setup:',t_langsetup,
166 & 'Stochastic MD step setup:',t_sdsetup,
168 write (iout,'(/28(1h=),a25,27(1h=))')
169 & ' End of MD calculation '
172 c-------------------------------------------------------------------------------
173 subroutine brown_step(itime)
174 c------------------------------------------------
175 c Perform a single Euler integration step of Brownian dynamics
176 c------------------------------------------------
177 implicit real*8 (a-h,o-z)
179 include 'COMMON.CONTROL'
182 include 'COMMON.LANGEVIN'
183 include 'COMMON.CHAIN'
184 include 'COMMON.DERIV'
186 include 'COMMON.LOCAL'
187 include 'COMMON.INTERACT'
188 include 'COMMON.IOUNITS'
189 include 'COMMON.NAMES'
190 include 'COMMON.TIME1'
191 double precision energia(0:n_ene),zapas(MAXRES6)
192 integer ilen,rstcount
194 double precision stochforcvec(MAXRES6)
195 double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
196 & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
197 & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
199 common /stochcalc/ stochforcvec
200 common /gucio/ cm, energia
202 logical lprn /.false./,lprn1 /.false./
204 double precision difftol /1.0d-5/
207 if (itype(i).ne.10) nbond=nbond+1
211 write (iout,*) "Generalized inverse of fricmat"
212 call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
224 Bmat(ind+j,ind1)=dC_norm(j,i)
229 if (itype(i).ne.10) then
232 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
238 write (iout,*) "Matrix Bmat"
239 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
245 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
250 write (iout,*) "Matrix GBmat"
251 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
257 Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
262 write (iout,*) "Matrix Cmat"
263 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
265 call matinvert(nbond,MAXRES2,Cmat,Cinv)
267 write (iout,*) "Matrix Cinv"
268 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
274 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
279 write (iout,*) "Matrix Tmat"
280 call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
290 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
295 write (iout,*) "Matrix Pmat"
296 call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
303 Td(i)=Td(i)+vbl*Tmat(i,ind)
306 if (itype(k).ne.10) then
308 Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind)
313 write (iout,*) "Vector Td"
315 write (iout,'(i5,f10.5)') i,Td(i)
318 call stochastic_force(stochforcvec)
320 write (iout,*) "stochforcvec"
322 write (iout,*) i,stochforcvec(i)
326 zapas(j)=-gcart(j,0)+stochforcvec(j)
328 dC_work(j)=dC_old(j,0)
334 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
335 dC_work(ind)=dC_old(j,i)
339 if (itype(i).ne.10) then
342 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
343 dC_work(ind)=dC_old(j,i+nres)
349 write (iout,*) "Initial d_t_work"
351 write (iout,*) i,d_t_work(i)
358 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
365 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
369 write (iout,*) "Final d_t_work and zapas"
371 write (iout,*) i,d_t_work(i),zapas(i)
385 dc_work(ind+j)=dc(j,i)
391 d_t(j,i+nres)=d_t_work(ind+j)
392 dc(j,i+nres)=zapas(ind+j)
393 dc_work(ind+j)=dc(j,i+nres)
399 write (iout,*) "Before correction for rotational lengthening"
400 write (iout,*) "New coordinates",
401 & " and differences between actual and standard bond lengths"
406 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
407 & i,(dC(j,i),j=1,3),xx
410 if (itype(i).ne.10) then
412 xx=vbld(i+nres)-vbldsc0(itype(i))
413 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
414 & i,(dC(j,i+nres),j=1,3),xx
418 c Second correction (rotational lengthening)
424 blen2 = scalar(dc(1,i),dc(1,i))
425 ppvec(ind)=2*vbl**2-blen2
426 diffbond=dabs(vbl-dsqrt(blen2))
427 if (diffbond.gt.diffmax) diffmax=diffbond
428 if (ppvec(ind).gt.0.0d0) then
429 ppvec(ind)=dsqrt(ppvec(ind))
434 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
438 if (itype(i).ne.10) then
440 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
441 ppvec(ind)=2*vbldsc0(itype(i))**2-blen2
442 diffbond=dabs(vbldsc0(itype(i))-dsqrt(blen2))
443 if (diffbond.gt.diffmax) diffmax=diffbond
444 if (ppvec(ind).gt.0.0d0) then
445 ppvec(ind)=dsqrt(ppvec(ind))
450 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
454 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
455 if (diffmax.lt.difftol) goto 10
459 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
465 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
476 dc_work(ind+j)=zapas(ind+j)
481 if (itype(i).ne.10) then
483 dc(j,i+nres)=zapas(ind+j)
484 dc_work(ind+j)=zapas(ind+j)
489 c Building the chain from the newly calculated coordinates
491 if (large.and. mod(itime,ntwe).eq.0) then
492 write (iout,*) "Cartesian and internal coordinates: step 1"
495 write (iout,'(a)') "Potential forces"
497 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
498 & (-gxcart(j,i),j=1,3)
500 write (iout,'(a)') "Stochastic forces"
502 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
503 & (stochforc(j,i+nres),j=1,3)
505 write (iout,'(a)') "Velocities"
507 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
508 & (d_t(j,i+nres),j=1,3)
512 write (iout,*) "After correction for rotational lengthening"
513 write (iout,*) "New coordinates",
514 & " and differences between actual and standard bond lengths"
519 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
520 & i,(dC(j,i),j=1,3),xx
523 if (itype(i).ne.10) then
525 xx=vbld(i+nres)-vbldsc0(itype(i))
526 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
527 & i,(dC(j,i+nres),j=1,3),xx
532 c write (iout,*) "Too many attempts at correcting the bonds"
536 c Calculate energy and forces
539 potE=energia(0)-energia(20)
542 c Calculate the kinetic and total energy and the kinetic temperature
544 t_enegrad=t_enegrad+tcpu()-tt0
546 kinetic_T=2.0d0/(dimen*Rb)*EK
549 c-------------------------------------------------------------------------------
550 subroutine velverlet_step(itime)
551 c-------------------------------------------------------------------------------
552 c Perform a single velocity Verlet step; the time step can be rescaled if
553 c increments in accelerations exceed the threshold
554 c-------------------------------------------------------------------------------
555 implicit real*8 (a-h,o-z)
557 include 'COMMON.CONTROL'
560 include 'COMMON.LANGEVIN'
561 include 'COMMON.CHAIN'
562 include 'COMMON.DERIV'
564 include 'COMMON.LOCAL'
565 include 'COMMON.INTERACT'
566 include 'COMMON.IOUNITS'
567 include 'COMMON.NAMES'
568 include 'COMMON.TIME1'
569 include 'COMMON.MUCA'
570 double precision energia(0:n_ene),vcm(3),incr(3)
571 double precision cm(3),L(3)
572 integer ilen,count,rstcount
575 integer maxcount_scale /20/
576 common /gucio/ cm, energia
577 double precision stochforcvec(MAXRES6)
578 common /stochcalc/ stochforcvec
586 else if (lang.eq.2 .or. lang.eq.3) then
587 call stochastic_force(stochforcvec)
591 icount_scale=icount_scale+1
592 if (icount_scale.gt.maxcount_scale) then
594 & "ERROR: too many attempts at scaling down the time step. ",
595 & "amax=",amax,"epdrift=",epdrift,
596 & "damax=",damax,"edriftmax=",edriftmax,
600 c First step of the velocity Verlet algorithm
603 else if (lang.eq.3) then
604 call sd_verlet1_ciccotti
605 else if (lang.eq.1) then
610 c Build the chain from the newly calculated coordinates
612 if (rattle) call rattle1
613 if (large.and. mod(itime,ntwe).eq.0) then
614 write (iout,*) "Cartesian and internal coordinates: step 1"
617 write (iout,*) "Accelerations"
619 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
620 & (d_a(j,i+nres),j=1,3)
622 write (iout,*) "Velocities, step 1"
624 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
625 & (d_t(j,i+nres),j=1,3)
629 c Calculate energy and forces
632 potE=energia(0)-energia(20)
634 c Get the new accelerations
636 t_enegrad=t_enegrad+tcpu()-tt0
637 c Determine maximum acceleration and scale down the timestep if needed
639 call predict_edrift(epdrift)
640 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
641 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
643 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
645 itime_scal=itime_scal+ifac_time
646 c fac_time=dmin1(damax/amax,0.5d0)
647 fac_time=0.5d0**ifac_time
648 d_time=d_time*fac_time
649 if (lang.eq.2 .or. lang.eq.3) then
650 c write (iout,*) "Calling sd_verlet_setup: 1"
651 c Rescale the stochastic forces and recalculate or restore
652 c the matrices of tinker integrator
653 if (itime_scal.gt.maxflag_stoch) then
654 if (large) write (iout,'(a,i5,a)')
655 & "Calculate matrices for stochastic step;",
656 & " itime_scal ",itime_scal
658 call sd_verlet_p_setup
660 call sd_verlet_ciccotti_setup
662 write (iout,'(2a,i3,a,i3,1h.)')
663 & "Warning: cannot store matrices for stochastic",
664 & " integration because the index",itime_scal,
665 & " is greater than",maxflag_stoch
666 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
667 & " integration Langevin algorithm for better efficiency."
668 else if (flag_stoch(itime_scal)) then
669 if (large) write (iout,'(a,i5,a,l1)')
670 & "Restore matrices for stochastic step; itime_scal ",
671 & itime_scal," flag ",flag_stoch(itime_scal)
674 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
675 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
676 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
677 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
678 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
679 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
683 if (large) write (iout,'(2a,i5,a,l1)')
684 & "Calculate & store matrices for stochastic step;",
685 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
687 call sd_verlet_p_setup
689 call sd_verlet_ciccotti_setup
691 flag_stoch(ifac_time)=.true.
694 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
695 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
696 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
697 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
698 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
699 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
703 fac_time=1.0d0/dsqrt(fac_time)
705 stochforcvec(i)=fac_time*stochforcvec(i)
707 else if (lang.eq.1) then
708 c Rescale the accelerations due to stochastic forces
709 fac_time=1.0d0/dsqrt(fac_time)
711 d_as_work(i)=d_as_work(i)*fac_time
714 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
715 & "itime",itime," Timestep scaled down to ",
716 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
718 c Second step of the velocity Verlet algorithm
721 else if (lang.eq.3) then
722 call sd_verlet2_ciccotti
723 else if (lang.eq.1) then
728 if (rattle) call rattle2
730 if (d_time.ne.d_time0) then
732 if (lang.eq.2 .or. lang.eq.3) then
733 if (large) write (iout,'(a)')
734 & "Restore original matrices for stochastic step"
735 c write (iout,*) "Calling sd_verlet_setup: 2"
736 c Restore the matrices of tinker integrator if the time step has been restored
739 pfric_mat(i,j)=pfric0_mat(i,j,0)
740 afric_mat(i,j)=afric0_mat(i,j,0)
741 vfric_mat(i,j)=vfric0_mat(i,j,0)
742 prand_mat(i,j)=prand0_mat(i,j,0)
743 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
744 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
752 c Calculate the kinetic and the total energy and the kinetic temperature
757 c write (iout,*) "step",itime," EK",EK," EK1",EK1
759 c Couple the system to Berendsen bath if needed
760 if (tbf .and. lang.eq.0) then
763 kinetic_T=2.0d0/(dimen*Rb)*EK
764 c Backup the coordinates, velocities, and accelerations
768 d_t_old(j,i)=d_t(j,i)
769 d_a_old(j,i)=d_a(j,i)
772 if (mod(itime,ntwe).eq.0 .and. large) then
773 write (iout,*) "Velocities, step 2"
775 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
776 & (d_t(j,i+nres),j=1,3)
781 c-------------------------------------------------------------------------------
782 subroutine RESPA_step(itime)
783 c-------------------------------------------------------------------------------
784 c Perform a single RESPA step.
785 c-------------------------------------------------------------------------------
786 implicit real*8 (a-h,o-z)
788 include 'COMMON.CONTROL'
791 include 'COMMON.LANGEVIN'
792 include 'COMMON.CHAIN'
793 include 'COMMON.DERIV'
795 include 'COMMON.LOCAL'
796 include 'COMMON.INTERACT'
797 include 'COMMON.IOUNITS'
798 include 'COMMON.NAMES'
799 include 'COMMON.TIME1'
800 double precision energia(0:n_ene),energia_short(0:n_ene),
801 & energia_long(0:n_ene)
802 double precision cm(3),L(3),vcm(3),incr(3)
803 integer ilen,count,rstcount
806 integer maxcount_scale /10/
808 double precision stochforcvec(MAXRES6)
809 common /stochcalc/ stochforcvec
812 if (large.and. mod(itime,ntwe).eq.0) then
813 write (iout,*) "***************** RESPA itime",itime
814 write (iout,*) "Cartesian and internal coordinates: step 0"
817 write (iout,*) "Accelerations"
819 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
820 & (d_a(j,i+nres),j=1,3)
822 write (iout,*) "Velocities, step 0"
824 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
825 & (d_t(j,i+nres),j=1,3)
829 c Perform the initial RESPA step (increment velocities)
830 c write (iout,*) "*********************** RESPA ini"
832 if (mod(itime,ntwe).eq.0 .and. large) then
833 write (iout,*) "Velocities, end"
835 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
836 & (d_t(j,i+nres),j=1,3)
839 c Compute the short-range forces
841 call etotal_short(energia_short)
847 d_t_old(j,i)=d_t(j,i)
848 d_a_old(j,i)=d_a(j,i)
852 c Split the time step
853 d_time=d_time/ntime_split
854 c Perform the short-range RESPSA steps (velocity Verlet increments of
855 c positions and velocities using short-range forces)
856 c write (iout,*) "*********************** RESPA split"
857 do itsplit=1,ntime_split
860 else if (lang.eq.2 .or. lang.eq.3) then
861 call stochastic_force(stochforcvec)
863 c First step of the velocity Verlet algorithm
866 else if (lang.eq.3) then
867 call sd_verlet1_ciccotti
868 else if (lang.eq.1) then
873 c Build the chain from the newly calculated coordinates
875 if (rattle) call rattle1
876 if (large.and. mod(itime,ntwe).eq.0) then
877 write (iout,*) "Cartesian and internal coordinates: step 1"
880 write (iout,*) "Accelerations"
882 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
883 & (d_a(j,i+nres),j=1,3)
885 write (iout,*) "Velocities, step 1"
887 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
888 & (d_t(j,i+nres),j=1,3)
892 c Calculate energy and forces
894 call etotal_short(energia_short)
896 c Get the new accelerations
898 t_enegrad=t_enegrad+tcpu()-tt0
899 c Second step of the velocity Verlet algorithm
902 else if (lang.eq.3) then
903 call sd_verlet2_ciccotti
904 else if (lang.eq.1) then
909 if (rattle) call rattle2
910 c Backup the coordinates, velocities, and accelerations
914 d_t_old(j,i)=d_t(j,i)
915 d_a_old(j,i)=d_a(j,i)
919 c Restore the time step
921 c Compute long-range forces
923 call etotal_long(energia_long)
925 c Compute accelerations from long-range forces
927 if (large.and. mod(itime,ntwe).eq.0) then
928 write (iout,*) "Cartesian and internal coordinates: step 2"
931 write (iout,*) "Accelerations"
933 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
934 & (d_a(j,i+nres),j=1,3)
936 write (iout,*) "Velocities, step 2"
938 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
939 & (d_t(j,i+nres),j=1,3)
942 c Compute the final RESPA step (increment velocities)
943 c write (iout,*) "*********************** RESPA fin"
945 c Compute the complete potential energy
946 potE=energia_short(0)+energia_long(0)
948 c Calculate the kinetic and the total energy and the kinetic temperature
951 c Couple the system to Berendsen bath if needed
952 if (tbf .and. lang.eq.0) then
955 kinetic_T=2.0d0/(dimen*Rb)*EK
956 c Backup the coordinates, velocities, and accelerations
957 if (mod(itime,ntwe).eq.0 .and. large) then
958 write (iout,*) "Velocities, end"
960 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
961 & (d_t(j,i+nres),j=1,3)
966 c---------------------------------------------------------------------
968 c First and last RESPA step (incrementing velocities using long-range
970 implicit real*8 (a-h,o-z)
972 include 'COMMON.CONTROL'
975 include 'COMMON.CHAIN'
976 include 'COMMON.DERIV'
978 include 'COMMON.LOCAL'
979 include 'COMMON.INTERACT'
980 include 'COMMON.IOUNITS'
981 include 'COMMON.NAMES'
983 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
987 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
991 if (itype(i).ne.10) then
994 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1000 c-----------------------------------------------------------------
1002 c Applying velocity Verlet algorithm - step 1 to coordinates
1003 implicit real*8 (a-h,o-z)
1004 include 'DIMENSIONS'
1005 include 'COMMON.CONTROL'
1006 include 'COMMON.VAR'
1008 include 'COMMON.CHAIN'
1009 include 'COMMON.DERIV'
1010 include 'COMMON.GEO'
1011 include 'COMMON.LOCAL'
1012 include 'COMMON.INTERACT'
1013 include 'COMMON.IOUNITS'
1014 include 'COMMON.NAMES'
1015 double precision adt,adt2
1018 adt=d_a_old(j,0)*d_time
1020 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1021 d_t_new(j,0)=d_t_old(j,0)+adt2
1022 d_t(j,0)=d_t_old(j,0)+adt
1026 adt=d_a_old(j,i)*d_time
1028 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1029 d_t_new(j,i)=d_t_old(j,i)+adt2
1030 d_t(j,i)=d_t_old(j,i)+adt
1034 if (itype(i).ne.10) then
1037 adt=d_a_old(j,inres)*d_time
1039 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1040 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1041 d_t(j,inres)=d_t_old(j,inres)+adt
1047 c---------------------------------------------------------------------
1049 c Step 2 of the velocity Verlet algorithm: update velocities
1050 implicit real*8 (a-h,o-z)
1051 include 'DIMENSIONS'
1052 include 'COMMON.CONTROL'
1053 include 'COMMON.VAR'
1055 include 'COMMON.CHAIN'
1056 include 'COMMON.DERIV'
1057 include 'COMMON.GEO'
1058 include 'COMMON.LOCAL'
1059 include 'COMMON.INTERACT'
1060 include 'COMMON.IOUNITS'
1061 include 'COMMON.NAMES'
1063 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1067 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1071 if (itype(i).ne.10) then
1074 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1080 c-----------------------------------------------------------------
1081 subroutine sddir_precalc
1082 c Applying velocity Verlet algorithm - step 1 to coordinates
1083 implicit real*8 (a-h,o-z)
1084 include 'DIMENSIONS'
1085 include 'COMMON.CONTROL'
1086 include 'COMMON.VAR'
1088 include 'COMMON.LANGEVIN'
1089 include 'COMMON.CHAIN'
1090 include 'COMMON.DERIV'
1091 include 'COMMON.GEO'
1092 include 'COMMON.LOCAL'
1093 include 'COMMON.INTERACT'
1094 include 'COMMON.IOUNITS'
1095 include 'COMMON.NAMES'
1096 double precision stochforcvec(MAXRES6)
1097 common /stochcalc/ stochforcvec
1099 c Compute friction and stochastic forces
1102 call stochastic_force(stochforcvec)
1104 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1105 c forces (d_as_work)
1111 d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1112 d_as_work(i)=d_as_work(i)+Ginv(i,j)*stochforcvec(j)
1117 c---------------------------------------------------------------------
1118 subroutine sddir_verlet1
1119 c Applying velocity Verlet algorithm - step 1 to velocities
1120 implicit real*8 (a-h,o-z)
1121 include 'DIMENSIONS'
1122 include 'COMMON.CONTROL'
1123 include 'COMMON.VAR'
1125 include 'COMMON.LANGEVIN'
1126 include 'COMMON.CHAIN'
1127 include 'COMMON.DERIV'
1128 include 'COMMON.GEO'
1129 include 'COMMON.LOCAL'
1130 include 'COMMON.INTERACT'
1131 include 'COMMON.IOUNITS'
1132 include 'COMMON.NAMES'
1133 c Revised 3/31/05 AL: correlation between random contributions to
1134 c position and velocity increments included.
1135 double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1136 double precision adt,adt2
1138 c Add the contribution from BOTH friction and stochastic force to the
1139 c coordinates, but ONLY the contribution from the friction forces to velocities
1142 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1143 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1144 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1145 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1146 d_t(j,0)=d_t_old(j,0)+adt
1151 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1152 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1153 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1154 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1155 d_t(j,i)=d_t_old(j,i)+adt
1160 if (itype(i).ne.10) then
1163 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1164 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1165 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1166 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1167 d_t(j,inres)=d_t_old(j,inres)+adt
1174 c---------------------------------------------------------------------
1175 subroutine sddir_verlet2
1176 c Calculating the adjusted velocities for accelerations
1177 implicit real*8 (a-h,o-z)
1178 include 'DIMENSIONS'
1179 include 'COMMON.CONTROL'
1180 include 'COMMON.VAR'
1182 include 'COMMON.LANGEVIN'
1183 include 'COMMON.CHAIN'
1184 include 'COMMON.DERIV'
1185 include 'COMMON.GEO'
1186 include 'COMMON.LOCAL'
1187 include 'COMMON.INTERACT'
1188 include 'COMMON.IOUNITS'
1189 include 'COMMON.NAMES'
1190 double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1191 double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1192 c Revised 3/31/05 AL: correlation between random contributions to
1193 c position and velocity increments included.
1194 c The correlation coefficients are calculated at low-friction limit.
1195 c Also, friction forces are now not calculated with new velocities.
1197 c call friction_force
1198 call stochastic_force(stochforcvec)
1200 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1201 c forces (d_as_work)
1204 c d_af_work(i)=0.0d0
1207 c d_af_work(i)=d_af_work(i)+Ginv(i,j)*fric_work(j)
1208 d_as_work1(i)=d_as_work1(i)+Ginv(i,j)*stochforcvec(j)
1215 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1216 & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1221 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1222 & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1227 if (itype(i).ne.10) then
1230 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1231 & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1232 & +cos60*d_as_work1(ind+j))*d_time
1239 c---------------------------------------------------------------------
1240 subroutine max_accel
1242 c Find the maximum difference in the accelerations of the the sites
1243 c at the beginning and the end of the time step.
1245 implicit real*8 (a-h,o-z)
1246 include 'DIMENSIONS'
1247 include 'COMMON.CONTROL'
1248 include 'COMMON.VAR'
1250 include 'COMMON.CHAIN'
1251 include 'COMMON.DERIV'
1252 include 'COMMON.GEO'
1253 include 'COMMON.LOCAL'
1254 include 'COMMON.INTERACT'
1255 include 'COMMON.IOUNITS'
1256 double precision aux(3),accel(3)
1258 aux(j)=d_a(j,0)-d_a_old(j,0)
1265 accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1266 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1273 if (itype(i).ne.10) then
1275 accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1279 if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1282 aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1287 c---------------------------------------------------------------------
1288 subroutine predict_edrift(epdrift)
1290 c Predict the drift of the potential energy
1292 implicit real*8 (a-h,o-z)
1293 include 'DIMENSIONS'
1294 include 'COMMON.CONTROL'
1295 include 'COMMON.VAR'
1297 include 'COMMON.CHAIN'
1298 include 'COMMON.DERIV'
1299 include 'COMMON.GEO'
1300 include 'COMMON.LOCAL'
1301 include 'COMMON.INTERACT'
1302 include 'COMMON.IOUNITS'
1303 include 'COMMON.MUCA'
1304 double precision epdrift,epdriftij
1305 c Drift of the potential energy
1311 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1312 if (lmuca) epdriftij=epdriftij*factor
1313 c write (iout,*) "back",i,j,epdriftij
1314 if (epdriftij.gt.epdrift) epdrift=epdriftij
1318 if (itype(i).ne.10) then
1321 & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1322 if (lmuca) epdriftij=epdriftij*factor
1323 c write (iout,*) "side",i,j,epdriftij
1324 if (epdriftij.gt.epdrift) epdrift=epdriftij
1328 epdrift=0.5d0*epdrift*d_time*d_time
1329 c write (iout,*) "epdrift",epdrift
1332 c-----------------------------------------------------------------------
1333 subroutine verlet_bath
1335 c Coupling to the thermostat by using the Berendsen algorithm
1337 implicit real*8 (a-h,o-z)
1338 include 'DIMENSIONS'
1339 include 'COMMON.CONTROL'
1340 include 'COMMON.VAR'
1342 include 'COMMON.CHAIN'
1343 include 'COMMON.DERIV'
1344 include 'COMMON.GEO'
1345 include 'COMMON.LOCAL'
1346 include 'COMMON.INTERACT'
1347 include 'COMMON.IOUNITS'
1348 include 'COMMON.NAMES'
1349 double precision T_half,fact
1351 T_half=2.0d0/(dimen*Rb)*EK
1352 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1353 c write(iout,*) "T_half", T_half
1354 c write(iout,*) "EK", EK
1355 c write(iout,*) "fact", fact
1357 d_t(j,0)=fact*d_t(j,0)
1361 d_t(j,i)=fact*d_t(j,i)
1365 if (itype(i).ne.10) then
1368 d_t(j,inres)=fact*d_t(j,inres)
1374 c---------------------------------------------------------
1376 c Set up the initial conditions of a MD simulation
1377 implicit real*8 (a-h,o-z)
1378 include 'DIMENSIONS'
1381 include 'COMMON.INFO'
1382 include 'COMMON.SETUP'
1386 include 'COMMON.CONTROL'
1387 include 'COMMON.VAR'
1389 include 'COMMON.LANGEVIN'
1390 include 'COMMON.CHAIN'
1391 include 'COMMON.DERIV'
1392 include 'COMMON.GEO'
1393 include 'COMMON.LOCAL'
1394 include 'COMMON.INTERACT'
1395 include 'COMMON.IOUNITS'
1396 include 'COMMON.NAMES'
1397 real*8 energia(0:n_ene),energia_long(0:n_ene),
1398 & energia_short(0:n_ene),vcm(3),incr(3)
1399 double precision cm(3),L(3),xv,sigv,lowb,highb
1406 write(iout,*) "d_time", d_time
1407 c Compute the standard deviations of stochastic forces for Langevin dynamics
1408 c if the friction coefficients do not depend on surface area
1409 if (lang.gt.0 .and. .not.surfarea) then
1411 stdforcp(i)=stdfp*dsqrt(gamp)
1414 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1417 c Open the pdb file for snapshotshots
1419 if (nprocs.eq.1) then
1422 npos = dlog10(dfloat(nprocs-1))+1
1424 if (npos.lt.3) npos=3
1425 write (liczba,'(i1)') npos
1426 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1428 write (liczba,form) myrank
1431 & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1434 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1439 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1441 cartname=prefix(:ilen(prefix))//"_MD.cx"
1445 write (qstr,'(256(1h ))')
1451 write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw
1452 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1460 write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw
1461 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1465 pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1466 cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1467 statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1471 write(iout,*) "Initial state will be read from file ",
1472 & rest2name(:ilen(rest2name))
1475 c Generate initial velocities
1476 write(iout,*) "Initial velocities randomly generated"
1480 c rest2name = prefix(:ilen(prefix))//'.rst'
1481 write(iout,*) "Initial backbone velocities"
1483 write(iout,*) (d_t(j,i),j=1,3)
1485 write(iout,*) "Initial side-chain velocities"
1487 write(iout,*) (d_t(j,i+nres),j=1,3)
1489 c Zeroing the total angular momentum of the system
1490 write(iout,*) "Calling the zero-angular
1491 & momentum subroutine"
1493 c Getting the potential energy and forces and velocities and accelerations
1495 c write (iout,*) "velocity of the center of the mass:"
1496 c write (iout,*) (vcm(j),j=1,3)
1498 d_t(j,0)=d_t(j,0)-vcm(j)
1500 c Removing the velocity of the center of mass
1502 write (iout,*) "vcm right after adjustment:"
1503 write (iout,*) (vcm(j),j=1,3)
1507 call chainbuild_cart
1510 call verlet_bath(EK)
1512 kinetic_T=2.0d0/(dimen*Rb)*EK
1516 call etotal(energia)
1521 if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1522 write(iout,*) "Potential energy"
1523 write(iout,*) (energia(i),i=0,n_ene)
1524 potE=energia(0)-energia(20)
1527 write (iout,*) "Initial:",
1528 & " Kinetic energy",EK," potential energy",potE,
1529 & " total energy",totE," maximum acceleration ",
1532 write (iout,*) "Initial velocities"
1534 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1535 & (d_t(j,i+nres),j=1,3)
1537 write (iout,*) "Initial accelerations"
1539 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1540 & (d_a(j,i+nres),j=1,3)
1546 d_t_old(j,i)=d_t(j,i)
1547 d_a_old(j,i)=d_a(j,i)
1549 c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1553 call etotal_long(energia_long)
1556 c call etotal_short(energia_short)
1557 c write (iout,*) "energia_long",energia_long(0),
1558 c & " energia_short",energia_short(0),
1559 c & " total",energia_long(0)+energia_short(0)
1563 c-----------------------------------------------------------
1564 subroutine random_vel
1565 implicit real*8 (a-h,o-z)
1566 include 'DIMENSIONS'
1567 include 'COMMON.CONTROL'
1568 include 'COMMON.VAR'
1570 include 'COMMON.LANGEVIN'
1571 include 'COMMON.CHAIN'
1572 include 'COMMON.DERIV'
1573 include 'COMMON.GEO'
1574 include 'COMMON.LOCAL'
1575 include 'COMMON.INTERACT'
1576 include 'COMMON.IOUNITS'
1577 include 'COMMON.NAMES'
1578 include 'COMMON.TIME1'
1579 double precision xv,sigv,lowb,highb
1580 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1581 c First generate velocities in the eigenspace of the G matrix
1582 c write (iout,*) "Calling random_vel"
1585 sigv=dsqrt((Rb*t_bath)/geigen(i))
1588 d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1592 c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1594 c Transform velocities to UNRES coordinate space
1598 d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1601 c Transfer to the d_t vector
1603 d_t(j,0)=d_t_work(j)
1609 d_t(j,i)=d_t_work(ind)
1613 if (itype(i).ne.10) then
1616 d_t(j,i+nres)=d_t_work(ind)
1621 c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1622 c & 2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1625 c-----------------------------------------------------------
1626 subroutine sd_verlet_p_setup
1627 c Sets up the parameters of stochastic Verlet algorithm
1628 implicit real*8 (a-h,o-z)
1629 include 'DIMENSIONS'
1630 include 'COMMON.CONTROL'
1631 include 'COMMON.VAR'
1633 include 'COMMON.LANGEVIN'
1634 include 'COMMON.CHAIN'
1635 include 'COMMON.DERIV'
1636 include 'COMMON.GEO'
1637 include 'COMMON.LOCAL'
1638 include 'COMMON.INTERACT'
1639 include 'COMMON.IOUNITS'
1640 include 'COMMON.NAMES'
1641 include 'COMMON.TIME1'
1642 double precision emgdt(MAXRES6),
1643 & pterm,vterm,rho,rhoc,vsig,
1644 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1645 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1646 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1647 logical lprn /.false./
1648 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1649 double precision ktm
1652 c AL 8/17/04 Code adapted from tinker
1654 c Get the frictional and random terms for stochastic dynamics in the
1655 c eigenspace of mass-scaled UNRES friction matrix
1658 gdt = fricgam(i) * d_time
1660 c Stochastic dynamics reduces to simple MD for zero friction
1662 if (gdt .le. zero) then
1663 pfric_vec(i) = 1.0d0
1664 vfric_vec(i) = d_time
1665 afric_vec(i) = 0.5d0 * d_time * d_time
1666 prand_vec(i) = 0.0d0
1667 vrand_vec1(i) = 0.0d0
1668 vrand_vec2(i) = 0.0d0
1670 c Analytical expressions when friction coefficient is large
1673 if (gdt .ge. gdt_radius) then
1676 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1677 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1678 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1679 vterm = 1.0d0 - egdt**2
1680 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1682 c Use series expansions when friction coefficient is small
1693 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1694 & - gdt5/120.0d0 + gdt6/720.0d0
1695 & - gdt7/5040.0d0 + gdt8/40320.0d0
1696 & - gdt9/362880.0d0) / fricgam(i)**2
1697 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1698 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1699 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1700 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1701 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1702 & + 127.0d0*gdt9/90720.0d0
1703 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1704 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1705 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1706 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1707 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1708 & - 17.0d0*gdt2/1280.0d0
1709 & + 17.0d0*gdt3/6144.0d0
1710 & + 40967.0d0*gdt4/34406400.0d0
1711 & - 57203.0d0*gdt5/275251200.0d0
1712 & - 1429487.0d0*gdt6/13212057600.0d0)
1715 c Compute the scaling factors of random terms for the nonzero friction case
1717 ktm = 0.5d0*d_time/fricgam(i)
1718 psig = dsqrt(ktm*pterm) / fricgam(i)
1719 vsig = dsqrt(ktm*vterm)
1720 rhoc = dsqrt(1.0d0 - rho*rho)
1722 vrand_vec1(i) = vsig * rho
1723 vrand_vec2(i) = vsig * rhoc
1728 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1731 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1732 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1736 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1738 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
1739 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
1740 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
1741 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
1742 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
1743 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
1744 c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
1745 c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
1746 c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
1747 c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
1748 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
1749 c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
1750 t_sdsetup=t_sdsetup+tcpu()-tt0
1753 c-------------------------------------------------------------
1754 subroutine eigtransf1(n,ndim,ab,d,c)
1757 double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1763 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1769 c-------------------------------------------------------------
1770 subroutine eigtransf(n,ndim,a,b,d,c)
1773 double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
1779 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
1785 c-------------------------------------------------------------
1786 subroutine sd_verlet1
1787 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
1788 implicit real*8 (a-h,o-z)
1789 include 'DIMENSIONS'
1790 include 'COMMON.CONTROL'
1791 include 'COMMON.VAR'
1793 include 'COMMON.LANGEVIN'
1794 include 'COMMON.CHAIN'
1795 include 'COMMON.DERIV'
1796 include 'COMMON.GEO'
1797 include 'COMMON.LOCAL'
1798 include 'COMMON.INTERACT'
1799 include 'COMMON.IOUNITS'
1800 include 'COMMON.NAMES'
1801 double precision stochforcvec(MAXRES6)
1802 common /stochcalc/ stochforcvec
1803 logical lprn /.false./
1805 c write (iout,*) "dc_old"
1807 c write (iout,'(i5,3f10.5,5x,3f10.5)')
1808 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
1811 dc_work(j)=dc_old(j,0)
1812 d_t_work(j)=d_t_old(j,0)
1813 d_a_work(j)=d_a_old(j,0)
1818 dc_work(ind+j)=dc_old(j,i)
1819 d_t_work(ind+j)=d_t_old(j,i)
1820 d_a_work(ind+j)=d_a_old(j,i)
1825 if (itype(i).ne.10) then
1827 dc_work(ind+j)=dc_old(j,i+nres)
1828 d_t_work(ind+j)=d_t_old(j,i+nres)
1829 d_a_work(ind+j)=d_a_old(j,i+nres)
1837 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
1841 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
1842 & vfric_mat(i,j),afric_mat(i,j),
1843 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
1851 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
1852 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
1853 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
1854 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
1856 d_t_work_new(i)=ddt1+0.5d0*ddt2
1857 d_t_work(i)=ddt1+ddt2
1861 d_t(j,0)=d_t_work(j)
1866 dc(j,i)=dc_work(ind+j)
1867 d_t(j,i)=d_t_work(ind+j)
1872 if (itype(i).ne.10) then
1875 dc(j,inres)=dc_work(ind+j)
1876 d_t(j,inres)=d_t_work(ind+j)
1883 c--------------------------------------------------------------------------
1884 subroutine sd_verlet2
1885 c Calculating the adjusted velocities for accelerations
1886 implicit real*8 (a-h,o-z)
1887 include 'DIMENSIONS'
1888 include 'COMMON.CONTROL'
1889 include 'COMMON.VAR'
1891 include 'COMMON.LANGEVIN'
1892 include 'COMMON.CHAIN'
1893 include 'COMMON.DERIV'
1894 include 'COMMON.GEO'
1895 include 'COMMON.LOCAL'
1896 include 'COMMON.INTERACT'
1897 include 'COMMON.IOUNITS'
1898 include 'COMMON.NAMES'
1899 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
1900 common /stochcalc/ stochforcvec
1902 c Compute the stochastic forces which contribute to velocity change
1904 call stochastic_force(stochforcvecV)
1910 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
1911 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
1912 & vrand_mat2(i,j)*stochforcvecV(j)
1914 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
1917 d_t(j,0)=d_t_work(j)
1922 d_t(j,i)=d_t_work(ind+j)
1927 if (itype(i).ne.10) then
1930 d_t(j,inres)=d_t_work(ind+j)
1937 c-----------------------------------------------------------
1938 subroutine sd_verlet_ciccotti_setup
1939 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
1941 implicit real*8 (a-h,o-z)
1942 include 'DIMENSIONS'
1943 include 'COMMON.CONTROL'
1944 include 'COMMON.VAR'
1946 include 'COMMON.LANGEVIN'
1947 include 'COMMON.CHAIN'
1948 include 'COMMON.DERIV'
1949 include 'COMMON.GEO'
1950 include 'COMMON.LOCAL'
1951 include 'COMMON.INTERACT'
1952 include 'COMMON.IOUNITS'
1953 include 'COMMON.NAMES'
1954 include 'COMMON.TIME1'
1955 double precision emgdt(MAXRES6),
1956 & pterm,vterm,rho,rhoc,vsig,
1957 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1958 & afric_vec(MAXRES6),prand_vec(MAXRES6),
1959 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1960 logical lprn /.false./
1961 double precision zero /1.0d-8/, gdt_radius /0.05d0/
1962 double precision ktm
1965 c AL 8/17/04 Code adapted from tinker
1967 c Get the frictional and random terms for stochastic dynamics in the
1968 c eigenspace of mass-scaled UNRES friction matrix
1971 write (iout,*) "i",i," fricgam",fricgam(i)
1972 gdt = fricgam(i) * d_time
1974 c Stochastic dynamics reduces to simple MD for zero friction
1976 if (gdt .le. zero) then
1977 pfric_vec(i) = 1.0d0
1978 vfric_vec(i) = d_time
1979 afric_vec(i) = 0.5d0*d_time*d_time
1980 prand_vec(i) = afric_vec(i)
1981 vrand_vec2(i) = vfric_vec(i)
1983 c Analytical expressions when friction coefficient is large
1988 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
1989 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
1990 prand_vec(i) = afric_vec(i)
1991 vrand_vec2(i) = vfric_vec(i)
1993 c Compute the scaling factors of random terms for the nonzero friction case
1995 c ktm = 0.5d0*d_time/fricgam(i)
1996 c psig = dsqrt(ktm*pterm) / fricgam(i)
1997 c vsig = dsqrt(ktm*vterm)
1998 c prand_vec(i) = psig*afric_vec(i)
1999 c vrand_vec2(i) = vsig*vfric_vec(i)
2004 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2007 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2008 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2012 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2014 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2015 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2016 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2017 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2018 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2019 t_sdsetup=t_sdsetup+tcpu()-tt0
2022 c-------------------------------------------------------------
2023 subroutine sd_verlet1_ciccotti
2024 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2025 implicit real*8 (a-h,o-z)
2026 include 'DIMENSIONS'
2027 include 'COMMON.CONTROL'
2028 include 'COMMON.VAR'
2030 include 'COMMON.LANGEVIN'
2031 include 'COMMON.CHAIN'
2032 include 'COMMON.DERIV'
2033 include 'COMMON.GEO'
2034 include 'COMMON.LOCAL'
2035 include 'COMMON.INTERACT'
2036 include 'COMMON.IOUNITS'
2037 include 'COMMON.NAMES'
2038 double precision stochforcvec(MAXRES6)
2039 common /stochcalc/ stochforcvec
2040 logical lprn /.false./
2042 c write (iout,*) "dc_old"
2044 c write (iout,'(i5,3f10.5,5x,3f10.5)')
2045 c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2048 dc_work(j)=dc_old(j,0)
2049 d_t_work(j)=d_t_old(j,0)
2050 d_a_work(j)=d_a_old(j,0)
2055 dc_work(ind+j)=dc_old(j,i)
2056 d_t_work(ind+j)=d_t_old(j,i)
2057 d_a_work(ind+j)=d_a_old(j,i)
2062 if (itype(i).ne.10) then
2064 dc_work(ind+j)=dc_old(j,i+nres)
2065 d_t_work(ind+j)=d_t_old(j,i+nres)
2066 d_a_work(ind+j)=d_a_old(j,i+nres)
2074 & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2078 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2079 & vfric_mat(i,j),afric_mat(i,j),
2080 & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2088 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2089 & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2090 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2091 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2093 d_t_work_new(i)=ddt1+0.5d0*ddt2
2094 d_t_work(i)=ddt1+ddt2
2098 d_t(j,0)=d_t_work(j)
2103 dc(j,i)=dc_work(ind+j)
2104 d_t(j,i)=d_t_work(ind+j)
2109 if (itype(i).ne.10) then
2112 dc(j,inres)=dc_work(ind+j)
2113 d_t(j,inres)=d_t_work(ind+j)
2120 c--------------------------------------------------------------------------
2121 subroutine sd_verlet2_ciccotti
2122 c Calculating the adjusted velocities for accelerations
2123 implicit real*8 (a-h,o-z)
2124 include 'DIMENSIONS'
2125 include 'COMMON.CONTROL'
2126 include 'COMMON.VAR'
2128 include 'COMMON.LANGEVIN'
2129 include 'COMMON.CHAIN'
2130 include 'COMMON.DERIV'
2131 include 'COMMON.GEO'
2132 include 'COMMON.LOCAL'
2133 include 'COMMON.INTERACT'
2134 include 'COMMON.IOUNITS'
2135 include 'COMMON.NAMES'
2136 double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2137 common /stochcalc/ stochforcvec
2139 c Compute the stochastic forces which contribute to velocity change
2141 call stochastic_force(stochforcvecV)
2147 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2148 c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2149 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2151 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2154 d_t(j,0)=d_t_work(j)
2159 d_t(j,i)=d_t_work(ind+j)
2164 if (itype(i).ne.10) then
2167 d_t(j,inres)=d_t_work(ind+j)