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

   1 c-------------------------------------------------------------------------------
   2       subroutine RESPA_step(itime)
   3 c-------------------------------------------------------------------------------
   4 c  Perform a single RESPA step.
   5 c-------------------------------------------------------------------------------
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.VAR'
  10       include 'COMMON.MD'
  11       include 'COMMON.LANGEVIN'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.DERIV'
  14       include 'COMMON.GEO'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.INTERACT'
  17       include 'COMMON.IOUNITS'
  18       include 'COMMON.NAMES'
  19       include 'COMMON.TIME1'
  20       double precision energia(0:n_ene),energia_short(0:n_ene),
  21      & energia_long(0:n_ene)
  22       double precision cm(3),L(3),vcm(3),incr(3)
  23       integer ilen,count,rstcount
  24       external ilen
  25       character*50 tytul
  26       integer maxcount_scale /10/
  27       common /gucio/ cm
  28       double precision stochforcvec(MAXRES6)
  29       double precision grad_tmp(3,0:maxres2)
  30       common /stochcalc/ stochforcvec
  31       integer itime
  32       logical scale
  33       double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
  34       if (large.and. mod(itime,ntwe).eq.0) then
  35         write (iout,*) "***************** RESPA itime",itime
  36         write (iout,*) "Cartesian and internal coordinates: step 0"
  37         call cartprint
  38         call intout
  39         write (iout,*) "Accelerations"
  40         do i=0,nres
  41           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
  42      &      (d_a(j,i+nres),j=1,3)
  43         enddo
  44         write (iout,*) "Velocities, step 0"
  45         do i=0,nres
  46           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  47      &      (d_t(j,i+nres),j=1,3)
  48         enddo
  49       endif
  50 c
  51 c Perform the initial RESPA step (increment velocities)
  52 c      write (iout,*) "*********************** RESPA ini"
  53
  54       if (tnp1) then
  55 creview          call tnp1_respa_step1
  56           call tnp_respa_step1
  57       else if (tnp) then
  58           call tnp_respa_step1
  59       else
  60           if (tnh.and..not.xiresp) then
  61             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
  62             do i=0,2*nres
  63              do j=1,3
  64               d_t(j,i)=d_t(j,i)*scale_nh
  65              enddo
  66             enddo
  67           endif
  68           call RESPA_vel
  69       endif
  70
  71 cd       if(tnp .or. tnp1) then
  72 cd        write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
  73 cd       endif
  74
  75       if (mod(itime,ntwe).eq.0 .and. large) then
  76         write (iout,*) "Velocities, end"
  77         do i=0,nres
  78           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  79      &      (d_t(j,i+nres),j=1,3)
  80         enddo
  81       endif
  82 c
  83 cr      if (tnp) then
  84 cr          do i=0,nres
  85 cr            do j=1,3
  86 cr             grad_tmp(j,i)=gcart(j,i)
  87 cr             grad_tmp(j,i+nres)=gxcart(j,i)
  88 cr            enddo
  89 cr          enddo
  90 cr      endif
  91 c
  92 c Compute the short-range forces
  93       call zerograd
  94       call etotal_short(energia_short)
  95       if (tnp.or.tnp1) potE=energia_short(0)
  96       call cartgrad
  97       call lagrangian
  98       do i=0,2*nres
  99         do j=1,3
 100           dc_old(j,i)=dc(j,i)
 101           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 102           d_a_old(j,i)=d_a(j,i)
 103         enddo
 104       enddo
 105       d_time0=d_time
 106 c Split the time step
 107       d_time=d_time/ntime_split
 108 ctest      E_long2=E_long
 109 c Perform the short-range RESPSA steps (velocity Verlet increments of
 110 c positions and velocities using short-range forces)
 111 c      write (iout,*) "*********************** RESPA split"
 112 creview      E_old=potE
 113       do itsplit=1,ntime_split
 114         if (lang.eq.1) then
 115           call sddir_precalc
 116         else if (lang.eq.2 .or. lang.eq.3) then
 117           call stochastic_force(stochforcvec)
 118         endif
 119 c First step of the velocity Verlet algorithm
 120         if (lang.eq.2) then
 121           call sd_verlet1
 122         else if (lang.eq.3) then
 123           call sd_verlet1_ciccotti
 124         else if (lang.eq.1) then
 125           call sddir_verlet1
 126         else if (tnp1) then
 127           call tnp1_respa_i_step1
 128 cr          if(itsplit.eq.1)then
 129 cr           d_time_s12=d_time0*0.5*s12_np
 130 cr           do j=1,3
 131 cr            d_t_half(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
 132 cr           enddo
 133 cr           do i=nnt,nct-1
 134 cr            do j=1,3
 135 cr             d_t_half(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
 136 cr            enddo
 137 cr           enddo
 138 cr           do i=nnt,nct
 139 cr            if (itype(i).ne.10) then
 140 cr             inres=i+nres
 141 cr             do j=1,3
 142 cr              d_t_half(j,inres)=d_t_old(j,inres)
 143 cr     &                +d_a_old(j,inres)*d_time_s12
 144 cr             enddo
 145 cr            endif
 146 cr           enddo
 147 cr          endif
 148         else if (tnp) then
 149           call tnp_respa_i_step1
 150         else
 151           if (tnh.and.xiresp) then
 152             call kinetic(EK)
 153             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
 154             do i=0,2*nres
 155              do j=1,3
 156               d_t_old(j,i)=d_t_old(j,i)*scale_nh
 157              enddo
 158             enddo
 159 cd            write(iout,*) "SSS1",itsplit,EK,scale_nh
 160           endif
 161           call verlet1
 162         endif
 163 c Build the chain from the newly calculated coordinates
 164         call chainbuild_cart
 165         if (rattle) call rattle1
 166         if (large.and. mod(itime,ntwe).eq.0) then
 167           write (iout,*) "Cartesian and internal coordinates: step 1"
 168           call cartprint
 169           call intout
 170           write (iout,*) "Accelerations"
 171           do i=0,nres
 172             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 173      &        (d_a(j,i+nres),j=1,3)
 174           enddo
 175           write (iout,*) "Velocities, step 1"
 176           do i=0,nres
 177             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 178      &        (d_t(j,i+nres),j=1,3)
 179           enddo
 180         endif
 181         tt0 = tcpu()
 182 c Calculate energy and forces
 183 c
 184 c E_long aproximation
 185 cr         if (tnp .or. tnp1) then
 186 cr          dtmp=0.5*d_time*(1.0/s12_np+1.0/s_np)
 187 cr          do i=0,2*nres
 188 cr            do j=1,3
 189 cr             E_long=E_long+d_t_new(j,i)*dtmp*grad_tmp(j,i)
 190 cr            enddo
 191 cr          enddo
 192 cr         endif
 193 c-------------------------------------
 194 c test of reviewer's comment
 195 cr        E_long=0
 196 c-------------------------------------
 197
 198
 199 c
 200 ctest        call etotal_long(energia_long)
 201 ctest        E_long=energia_long(0)
 202 ctest
 203         call zerograd
 204
 205         call etotal_short(energia_short)
 206         E_old=potE
 207         potE=energia_short(0)
 208
 209 c       if(tnp .or. tnp1) then
 210 c        write (iout,*) "kkk",E_long2,E_long
 211 c       endif
 212
 213
 214         call cartgrad
 215 c Get the new accelerations
 216         call lagrangian
 217         t_enegrad=t_enegrad+tcpu()-tt0
 218 c Second step of the velocity Verlet algorithm
 219         if (lang.eq.2) then
 220           call sd_verlet2
 221         else if (lang.eq.3) then
 222           call sd_verlet2_ciccotti
 223         else if (lang.eq.1) then
 224           call sddir_verlet2
 225         else if (tnp1) then
 226             call tnp1_respa_i_step2
 227         else if (tnp) then
 228             call tnp_respa_i_step2
 229         else
 230           call verlet2
 231           if (tnh.and.xiresp) then
 232             call kinetic(EK)
 233             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
 234             do i=0,2*nres
 235              do j=1,3
 236               d_t(j,i)=d_t(j,i)*scale_nh
 237              enddo
 238             enddo
 239 cd            write(iout,*) "SSS2",itsplit,EK,scale_nh
 240           endif
 241
 242         endif
 243         if (rattle) call rattle2
 244 c Backup the coordinates, velocities, and accelerations
 245         if (tnp .or. tnp1) then
 246          do i=0,2*nres
 247           do j=1,3
 248             d_t_old(j,i)=d_t(j,i)
 249             if (tnp) d_t(j,i)=d_t(j,i)/s_np
 250             if (tnp1) d_t(j,i)=d_t(j,i)/s_np
 251           enddo
 252          enddo
 253
 254         endif
 255         do i=0,2*nres
 256           do j=1,3
 257             dc_old(j,i)=dc(j,i)
 258             if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 259             d_a_old(j,i)=d_a(j,i)
 260           enddo
 261         enddo
 262
 263 cd       if(tnp .or. tnp1) then
 264 cd        call kinetic(EK)
 265 cd        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
 266 cd        H=(HNose1-H0)*s_np
 267 cd        write (iout,*) "jjj",EK,potE
 268 cd        write (iout,*) "iii H=",H,abs(HNose1-H0)/H0
 269 cd        write (iout,'(a,3f)')
 270 cd     &             "III NP S, pi",totT+itsplit*d_time, s_np, pi_np
 271 cd       endif
 272
 273
 274       enddo
 275 c Restore the time step
 276       d_time=d_time0
 277 c Compute long-range forces
 278       call zerograd
 279       call etotal_long(energia_long)
 280       E_long=energia_long(0)
 281       potE=energia_short(0)+energia_long(0)
 282       call cartgrad
 283 c Compute accelerations from long-range forces
 284       call lagrangian
 285       if (large.and. mod(itime,ntwe).eq.0) then
 286         write (iout,*) "Cartesian and internal coordinates: step 2"
 287         call cartprint
 288         call intout
 289         write (iout,*) "Accelerations"
 290         do i=0,nres
 291           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 292      &      (d_a(j,i+nres),j=1,3)
 293         enddo
 294         write (iout,*) "Velocities, step 2"
 295         do i=0,nres
 296           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 297      &      (d_t(j,i+nres),j=1,3)
 298         enddo
 299       endif
 300 c Compute the final RESPA step (increment velocities)
 301 c      write (iout,*) "*********************** RESPA fin"
 302       if (tnp1) then
 303 creview          call tnp1_respa_step2
 304           call tnp_respa_step2
 305       else if (tnp) then
 306           call tnp_respa_step2
 307       else
 308           call RESPA_vel
 309           if (tnh.and..not.xiresp) then
 310             call kinetic(EK)
 311             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
 312             do i=0,2*nres
 313              do j=1,3
 314               d_t(j,i)=d_t(j,i)*scale_nh
 315              enddo
 316             enddo
 317           endif
 318       endif
 319
 320         if (tnp .or. tnp1) then
 321          do i=0,2*nres
 322           do j=1,3
 323             d_t(j,i)=d_t_old(j,i)/s_np
 324           enddo
 325          enddo
 326         endif
 327
 328 c Compute the complete potential energy
 329 cc      potE=energia_short(0)+energia_long(0)
 330       totT=totT+d_time
 331 c Calculate the kinetic and the total energy and the kinetic temperature
 332       call kinetic(EK)
 333       totE=EK+potE
 334
 335 c Couple the system to Berendsen bath if needed
 336       if (tbf .and. lang.eq.0) then
 337         call verlet_bath
 338       endif
 339       kinetic_T=2.0d0/(dimen*Rb)*EK
 340 c Backup the coordinates, velocities, and accelerations
 341       if (mod(itime,ntwe).eq.0 .and. large) then
 342         write (iout,*) "Velocities, end"
 343         do i=0,nres
 344           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 345      &      (d_t(j,i+nres),j=1,3)
 346         enddo
 347       endif
 348
 349 c-----review
 350 c       if(tnp .or. tnp1) then
 351 c        HNose1=Hnose(EK,s_np,energia_short(0),pi_np,Q_np,t_bath,dimen)
 352 c_new_var_csplit         Csplit=H0-E_long
 353 c         Csplit=H0-energia_short(0)
 354 c       endif
 355 c----------
 356
 357
 358       if (mod(itime,ntwe).eq.0) then
 359
 360        if(tnp .or. tnp1) then
 361         write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
 362      &                          E_long,energia_short(0)
 363         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
 364         H=(HNose1-H0)*s_np
 365 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
 366 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
 367         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 368 cd        write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
 369        endif
 370
 371        if(tnh) then
 372         HNose1=Hnose_nh(EK,potE)
 373         H=HNose1-H0
 374         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 375        endif
 376
 377
 378        if (large) then
 379        itnp=0
 380        do j=1,3
 381         itnp=itnp+1
 382         vtnp(itnp)=d_t(j,0)
 383        enddo
 384        do i=nnt,nct-1
 385         do j=1,3
 386           itnp=itnp+1
 387           vtnp(itnp)=d_t(j,i)
 388         enddo
 389        enddo
 390        do i=nnt,nct
 391         if (itype(i).ne.10) then
 392           inres=i+nres
 393           do j=1,3
 394            itnp=itnp+1
 395            vtnp(itnp)=d_t(j,inres)
 396           enddo
 397         endif
 398        enddo
 399
 400 c Transform velocities from UNRES coordinate space to cartesian and Gvec
 401 c eigenvector space
 402
 403         do i=1,dimen
 404           vtnp_(i)=0.0d0
 405           vtnp_a(i)=0.0d0
 406           do j=1,dimen
 407             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
 408             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
 409           enddo
 410           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
 411         enddo
 412
 413         do i=1,dimen
 414          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
 415         enddo
 416
 417        endif
 418       endif
 419       return
 420       end