+++ /dev/null
- subroutine MD
-c------------------------------------------------
-c The driver for molecular dynamics subroutines
-c------------------------------------------------
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
- include 'COMMON.LANGEVIN'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision cm(3),L(3),vcm(3)
- double precision energia(0:n_ene)
- integer ilen,rstcount
- external ilen
- character*50 tytul
- common /gucio/ cm,energia
- integer itime
-c
- t_MDsetup=0.0d0
- t_langsetup=0.0d0
- t_MD=0.0d0
- t_enegrad=0.0d0
- t_sdsetup=0.0d0
- write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
- tt0 = tcpu()
-c Determine the inverse of the inertia matrix.
- call setup_MD_matrices
-c Initialize MD
- call init_MD
- t_MDsetup = tcpu()-tt0
- rstcount=0
-c Entering the MD loop
- tt0 = tcpu()
- if (lang.eq.2 .or. lang.eq.3) then
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
- do i=1,dimen
- do j=1,dimen
- pfric0_mat(i,j,0)=pfric_mat(i,j)
- afric0_mat(i,j,0)=afric_mat(i,j)
- vfric0_mat(i,j,0)=vfric_mat(i,j)
- prand0_mat(i,j,0)=prand_mat(i,j)
- vrand0_mat1(i,j,0)=vrand_mat1(i,j)
- vrand0_mat2(i,j,0)=vrand_mat2(i,j)
- enddo
- enddo
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
- else if (lang.eq.1 .or. lang.eq.4) then
- call setup_fricmat
- endif
- t_langsetup=tcpu()-tt0
- tt0=tcpu()
- do itime=1,n_timestep
- rstcount=rstcount+1
- if (lang.gt.0 .and. surfarea .and.
- & mod(itime,reset_fricmat).eq.0) then
- if (lang.eq.2 .or. lang.eq.3) then
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
- do i=1,dimen
- do j=1,dimen
- pfric0_mat(i,j,0)=pfric_mat(i,j)
- afric0_mat(i,j,0)=afric_mat(i,j)
- vfric0_mat(i,j,0)=vfric_mat(i,j)
- prand0_mat(i,j,0)=prand_mat(i,j)
- vrand0_mat1(i,j,0)=vrand_mat1(i,j)
- vrand0_mat2(i,j,0)=vrand_mat2(i,j)
- enddo
- enddo
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
- else if (lang.eq.1 .or. lang.eq.4) then
- call setup_fricmat
- endif
- write (iout,'(a,i10)')
- & "Friction matrix reset based on surface area, itime",itime
- endif
- if (reset_vel .and. tbf .and. lang.eq.0
- & .and. mod(itime,count_reset_vel).eq.0) then
- call random_vel
- write(iout,'(a,f20.2)')
- & "Velocities reset to random values, time",totT
- do i=0,2*nres
- do j=1,3
- d_t_old(j,i)=d_t(j,i)
- enddo
- enddo
- endif
- if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
- call inertia_tensor
- call vcm_vel(vcm)
- do j=1,3
- d_t(j,0)=d_t(j,0)-vcm(j)
- enddo
- call kinetic(EK)
- kinetic_T=2.0d0/(dimen*Rb)*EK
- scalfac=dsqrt(T_bath/kinetic_T)
- write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
- do i=0,2*nres
- do j=1,3
- d_t_old(j,i)=scalfac*d_t(j,i)
- enddo
- enddo
- endif
- if (lang.ne.4) then
- if (RESPA) then
-c Time-reversible RESPA algorithm
-c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
- call RESPA_step(itime)
- else
-c Variable time step algorithm.
- call velverlet_step(itime)
- endif
- else
- call brown_step(itime)
- endif
- if (mod(itime,ntwe).eq.0) call statout(itime)
- if (mod(itime,ntwx).eq.0) then
- write (tytul,'("time",f8.2)') totT
- if(mdpdb) then
- call pdbout(potE,tytul,ipdb)
- else
- call cartout(totT)
- endif
- endif
- if (rstcount.eq.1000.or.itime.eq.n_timestep) then
- open(irest2,file=rest2name,status='unknown')
- write(irest2,*) totT,EK,potE,totE
- do i=1,2*nres
- write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
- enddo
- do i=1,2*nres
- write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
- enddo
- close(irest2)
- rstcount=0
- endif
- enddo
- t_MD=tcpu()-tt0
- write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
- & ' Timing ',
- & 'MD calculations setup:',t_MDsetup,
- & 'Energy & gradient evaluation:',t_enegrad,
- & 'Stochastic MD setup:',t_langsetup,
- & 'Stochastic MD step setup:',t_sdsetup,
- & 'MD steps:',t_MD
- write (iout,'(/28(1h=),a25,27(1h=))')
- & ' End of MD calculation '
- return
- end
projects / unres.git / blobdiff