source/unres/src_MD/md-diff/np/tnp_step1.f

   1 c---------------------------------------------------------------------
   2       subroutine tnp_step1
   3 c Applying Nose-Poincare algorithm - step 1 to coordinates
   4 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
   5 c
   6 c d_t is not updated here, it is destroyed
   7 c
   8       implicit real*8 (a-h,o-z)
   9       include 'DIMENSIONS'
  10       include 'COMMON.CONTROL'
  11       include 'COMMON.VAR'
  12       include 'COMMON.MD'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.DERIV'
  15       include 'COMMON.GEO'
  16       include 'COMMON.LOCAL'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.NAMES'
  20       double precision C_np,d_time_s,tmp,d_time_ss
  21
  22       d_time_s=d_time*0.5*s_np
  23
  24       do j=1,3
  25         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
  26       enddo
  27       do i=nnt,nct-1
  28         do j=1,3
  29           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
  30         enddo
  31       enddo
  32       do i=nnt,nct
  33         if (itype(i).ne.10) then
  34           inres=i+nres
  35           do j=1,3
  36            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
  37           enddo
  38         endif
  39       enddo
  40
  41       do i=0,2*nres
  42        do j=1,3
  43         d_t(j,i)=d_t_new(j,i)
  44        enddo
  45       enddo
  46
  47       call kinetic(EK)
  48       EK=EK/s_np**2
  49
  50       C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
  51      &                     -pi_np
  52
  53       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
  54       tmp=0.5*d_time*pistar/Q_np
  55       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
  56 c      write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
  57
  58       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
  59
  60       do j=1,3
  61         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
  62       enddo
  63       do i=nnt,nct-1
  64         do j=1,3
  65           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
  66         enddo
  67       enddo
  68       do i=nnt,nct
  69         if (itype(i).ne.10) then
  70           inres=i+nres
  71           do j=1,3
  72            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
  73           enddo
  74         endif
  75       enddo
  76
  77       return
  78       end