source/unres/src_MD/md-diff/np/md.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       include 'COMMON.CONTROL'
   8       include 'COMMON.VAR'
   9       include 'COMMON.MD'
  10       include 'COMMON.LANGEVIN'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.DERIV'
  13       include 'COMMON.GEO'
  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)
  21       integer ilen,rstcount
  22       external ilen
  23       character*50 tytul
  24       common /gucio/ cm,energia
  25       integer itime
  26 c
  27       t_MDsetup=0.0d0
  28       t_langsetup=0.0d0
  29       t_MD=0.0d0
  30       t_enegrad=0.0d0
  31       t_sdsetup=0.0d0
  32       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  33       tt0 = tcpu()
  34 c Determine the inverse of the inertia matrix.
  35       call setup_MD_matrices
  36 c Initialize MD
  37       call init_MD
  38       t_MDsetup = tcpu()-tt0
  39       rstcount=0
  40 c   Entering the MD loop
  41       tt0 = tcpu()
  42       if (lang.eq.2 .or. lang.eq.3) then
  43         call setup_fricmat
  44         if (lang.eq.2) then
  45           call sd_verlet_p_setup
  46         else
  47           call sd_verlet_ciccotti_setup
  48         endif
  49         do i=1,dimen
  50           do j=1,dimen
  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)
  57           enddo
  58         enddo
  59         flag_stoch(0)=.true.
  60         do i=1,maxflag_stoch
  61           flag_stoch(i)=.false.
  62         enddo
  63       else if (lang.eq.1 .or. lang.eq.4) then
  64         call setup_fricmat
  65       endif
  66       t_langsetup=tcpu()-tt0
  67       tt0=tcpu()
  68       do itime=1,n_timestep
  69         rstcount=rstcount+1
  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
  73             call setup_fricmat
  74             if (lang.eq.2) then
  75               call sd_verlet_p_setup
  76             else
  77               call sd_verlet_ciccotti_setup
  78             endif
  79             do i=1,dimen
  80               do j=1,dimen
  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)
  87               enddo
  88             enddo
  89             flag_stoch(0)=.true.
  90             do i=1,maxflag_stoch
  91               flag_stoch(i)=.false.
  92             enddo
  93           else if (lang.eq.1 .or. lang.eq.4) then
  94             call setup_fricmat
  95           endif
  96           write (iout,'(a,i10)')
  97      &      "Friction matrix reset based on surface area, itime",itime
  98         endif
  99         if (reset_vel .and. tbf .and. lang.eq.0
 100      &      .and. mod(itime,count_reset_vel).eq.0) then
 101           call random_vel
 102           write(iout,'(a,f20.2)')
 103      &     "Velocities reset to random values, time",totT
 104           do i=0,2*nres
 105             do j=1,3
 106               d_t_old(j,i)=d_t(j,i)
 107             enddo
 108           enddo
 109         endif
 110         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 111           call inertia_tensor
 112           call vcm_vel(vcm)
 113           do j=1,3
 114              d_t(j,0)=d_t(j,0)-vcm(j)
 115           enddo
 116           call kinetic(EK)
 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
 120           do i=0,2*nres
 121             do j=1,3
 122               d_t_old(j,i)=scalfac*d_t(j,i)
 123             enddo
 124           enddo
 125         endif
 126         if (lang.ne.4) then
 127           if (RESPA) then
 128 c Time-reversible RESPA algorithm
 129 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 130             call RESPA_step(itime)
 131           else
 132 c Variable time step algorithm.
 133             call velverlet_step(itime)
 134           endif
 135         else
 136           call brown_step(itime)
 137         endif
 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
 141           if(mdpdb) then
 142              call pdbout(potE,tytul,ipdb)
 143           else
 144              call cartout(totT)
 145           endif
 146         endif
 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
 150            do i=1,2*nres
 151             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 152            enddo
 153            do i=1,2*nres
 154             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 155            enddo
 156           close(irest2)
 157           rstcount=0
 158         endif
 159       enddo
 160       t_MD=tcpu()-tt0
 161       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 162      &  '  Timing  ',
 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,
 167      & 'MD steps:',t_MD
 168       write (iout,'(/28(1h=),a25,27(1h=))')
 169      & '  End of MD calculation  '
 170       return
 171       end