source/unres/src_MD/MD_A-MTS.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 #ifdef MPI
   8       include "mpif.h"
   9       integer IERROR,ERRCODE
  10 #endif
  11       include 'COMMON.SETUP'
  12       include 'COMMON.CONTROL'
  13       include 'COMMON.VAR'
  14       include 'COMMON.MD'
  15 #ifndef LANG0
  16       include 'COMMON.LANGEVIN'
  17 #else
  18       include 'COMMON.LANGEVIN.lang0'
  19 #endif
  20       include 'COMMON.CHAIN'
  21       include 'COMMON.DERIV'
  22       include 'COMMON.GEO'
  23       include 'COMMON.LOCAL'
  24       include 'COMMON.INTERACT'
  25       include 'COMMON.IOUNITS'
  26       include 'COMMON.NAMES'
  27       include 'COMMON.TIME1'
  28       double precision cm(3),L(3),vcm(3)
  29 #ifdef VOUT
  30       double precision v_work(maxres6),v_transf(maxres6)
  31 #endif
  32       integer ilen,rstcount
  33       external ilen
  34       character*50 tytul
  35       common /gucio/ cm
  36       integer itime
  37 c
  38 #ifdef MPI
  39       if (ilen(tmpdir).gt.0)
  40      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
  41      &        //liczba(:ilen(liczba))//'.rst')
  42 #else
  43       if (ilen(tmpdir).gt.0)
  44      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
  45 #endif
  46       t_MDsetup=0.0d0
  47       t_langsetup=0.0d0
  48       t_MD=0.0d0
  49       t_enegrad=0.0d0
  50       t_sdsetup=0.0d0
  51       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  52 #ifdef MPI
  53       tt0=MPI_Wtime()
  54 #else
  55       tt0 = tcpu()
  56 #endif
  57 c Determine the inverse of the inertia matrix.
  58       call setup_MD_matrices
  59 c Initialize MD
  60       call init_MD
  61 #ifdef MPI
  62       t_MDsetup = MPI_Wtime()-tt0
  63 #else
  64       t_MDsetup = tcpu()-tt0
  65 #endif
  66       rstcount=0
  67 c   Entering the MD loop
  68 #ifdef MPI
  69       tt0 = MPI_Wtime()
  70 #else
  71       tt0 = tcpu()
  72 #endif
  73       if (lang.eq.2 .or. lang.eq.3) then
  74 #ifndef   LANG0
  75         call setup_fricmat
  76         if (lang.eq.2) then
  77           call sd_verlet_p_setup
  78         else
  79           call sd_verlet_ciccotti_setup
  80         endif
  81         do i=1,dimen3
  82           do j=1,dimen3
  83             pfric0_mat(i,j,0)=pfric_mat(i,j)
  84             afric0_mat(i,j,0)=afric_mat(i,j)
  85             vfric0_mat(i,j,0)=vfric_mat(i,j)
  86             prand0_mat(i,j,0)=prand_mat(i,j)
  87             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
  88             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
  89           enddo
  90         enddo
  91         flag_stoch(0)=.true.
  92         do i=1,maxflag_stoch
  93           flag_stoch(i)=.false.
  94         enddo
  95 #else
  96         write (iout,*)
  97      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
  98 #ifdef MPI
  99         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 100 #endif
 101         stop
 102 #endif
 103       else if (lang.eq.1 .or. lang.eq.4) then
 104         call setup_fricmat
 105       endif
 106 #ifdef MPI
 107       t_langsetup=MPI_Wtime()-tt0
 108       tt0=MPI_Wtime()
 109 #else
 110       t_langsetup=tcpu()-tt0
 111       tt0=tcpu()
 112 #endif
 113       do itime=1,n_timestep
 114         rstcount=rstcount+1
 115         if (lang.gt.0 .and. surfarea .and.
 116      &      mod(itime,reset_fricmat).eq.0) then
 117           if (lang.eq.2 .or. lang.eq.3) then
 118 #ifndef LANG0
 119             call setup_fricmat
 120             if (lang.eq.2) then
 121               call sd_verlet_p_setup
 122             else
 123               call sd_verlet_ciccotti_setup
 124             endif
 125             do i=1,dimen3
 126               do j=1,dimen3
 127                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 128                 afric0_mat(i,j,0)=afric_mat(i,j)
 129                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 130                 prand0_mat(i,j,0)=prand_mat(i,j)
 131                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 132                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 133               enddo
 134             enddo
 135             flag_stoch(0)=.true.
 136             do i=1,maxflag_stoch
 137               flag_stoch(i)=.false.
 138             enddo
 139 #endif
 140           else if (lang.eq.1 .or. lang.eq.4) then
 141             call setup_fricmat
 142           endif
 143           write (iout,'(a,i10)')
 144      &      "Friction matrix reset based on surface area, itime",itime
 145         endif
 146         if (reset_vel .and. tbf .and. lang.eq.0
 147      &      .and. mod(itime,count_reset_vel).eq.0) then
 148           call random_vel
 149           write(iout,'(a,f20.2)')
 150      &     "Velocities reset to random values, time",totT
 151           do i=0,2*nres
 152             do j=1,3
 153               d_t_old(j,i)=d_t(j,i)
 154             enddo
 155           enddo
 156         endif
 157         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 158           call inertia_tensor
 159           call vcm_vel(vcm)
 160           do j=1,3
 161              d_t(j,0)=d_t(j,0)-vcm(j)
 162           enddo
 163           call kinetic(EK)
 164           kinetic_T=2.0d0/(dimen3*Rb)*EK
 165           scalfac=dsqrt(T_bath/kinetic_T)
 166           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 167           do i=0,2*nres
 168             do j=1,3
 169               d_t_old(j,i)=scalfac*d_t(j,i)
 170             enddo
 171           enddo
 172         endif
 173         if (lang.ne.4) then
 174           if (RESPA) then
 175 c Time-reversible RESPA algorithm
 176 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 177             call RESPA_step(itime)
 178           else
 179 c Variable time step algorithm.
 180             call velverlet_step(itime)
 181           endif
 182         else
 183 #ifdef BROWN
 184           call brown_step(itime)
 185 #else
 186           print *,"Brown dynamics not here!"
 187 #ifdef MPI
 188           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 189 #endif
 190           stop
 191 #endif
 192         endif
 193         if (ntwe.ne.0) then
 194          if (mod(itime,ntwe).eq.0) call statout(itime)
 195 #ifdef VOUT
 196         do j=1,3
 197           v_work(j)=d_t(j,0)
 198         enddo
 199         ind=3
 200         do i=nnt,nct-1
 201           do j=1,3
 202             ind=ind+1
 203             v_work(ind)=d_t(j,i)
 204           enddo
 205         enddo
 206         do i=nnt,nct
 207           if (itype(i).ne.10) then
 208             do j=1,3
 209               ind=ind+1
 210               v_work(ind)=d_t(j,i+nres)
 211             enddo
 212           endif
 213         enddo
 214
 215         write (66,'(80f10.5)')
 216      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 217         do i=1,ind
 218           v_transf(i)=0.0d0
 219           do j=1,ind
 220             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 221           enddo
 222            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 223         enddo
 224         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 225 #endif
 226         endif
 227         if (mod(itime,ntwx).eq.0) then
 228           write (tytul,'("time",f8.2)') totT
 229           if(mdpdb) then
 230              call pdbout(potE,tytul,ipdb)
 231           else
 232              call cartout(totT)
 233           endif
 234         endif
 235         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
 236            open(irest2,file=rest2name,status='unknown')
 237            write(irest2,*) totT,EK,potE,totE,t_bath
 238            do i=1,2*nres
 239             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 240            enddo
 241            do i=1,2*nres
 242             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 243            enddo
 244           close(irest2)
 245           rstcount=0
 246         endif
 247       enddo
 248 #ifdef MPI
 249       t_MD=MPI_Wtime()-tt0
 250 #else
 251       t_MD=tcpu()-tt0
 252 #endif
 253       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 254      &  '  Timing  ',
 255      & 'MD calculations setup:',t_MDsetup,
 256      & 'Energy & gradient evaluation:',t_enegrad,
 257      & 'Stochastic MD setup:',t_langsetup,
 258      & 'Stochastic MD step setup:',t_sdsetup,
 259      & 'MD steps:',t_MD
 260       write (iout,'(/28(1h=),a25,27(1h=))')
 261      & '  End of MD calculation  '
 262 #ifdef TIMING_ENE
 263       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
 264      &  " eshort",t_eshort
 265       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
 266      & " time_fricmatmult",time_fricmatmult," time_fsample ",
 267      & time_fsample
 268 #endif
 269       return
 270       end
 271 c-------------------------------------------------------------------------------
 272       subroutine velverlet_step(itime)
 273 c-------------------------------------------------------------------------------
 274 c  Perform a single velocity Verlet step; the time step can be rescaled if
 275 c  increments in accelerations exceed the threshold
 276 c-------------------------------------------------------------------------------
 277       implicit real*8 (a-h,o-z)
 278       include 'DIMENSIONS'
 279 #ifdef MPI
 280       include 'mpif.h'
 281       integer ierror,ierrcode
 282 #endif
 283       include 'COMMON.SETUP'
 284       include 'COMMON.CONTROL'
 285       include 'COMMON.VAR'
 286       include 'COMMON.MD'
 287 #ifndef LANG0
 288       include 'COMMON.LANGEVIN'
 289 #else
 290       include 'COMMON.LANGEVIN.lang0'
 291 #endif
 292       include 'COMMON.CHAIN'
 293       include 'COMMON.DERIV'
 294       include 'COMMON.GEO'
 295       include 'COMMON.LOCAL'
 296       include 'COMMON.INTERACT'
 297       include 'COMMON.IOUNITS'
 298       include 'COMMON.NAMES'
 299       include 'COMMON.TIME1'
 300       include 'COMMON.MUCA'
 301       double precision vcm(3),incr(3)
 302       double precision cm(3),L(3)
 303       integer ilen,count,rstcount
 304       external ilen
 305       character*50 tytul
 306       integer maxcount_scale /20/
 307       common /gucio/ cm
 308       double precision stochforcvec(MAXRES6)
 309       common /stochcalc/ stochforcvec
 310       integer itime
 311       logical scale
 312       double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
 313       double precision vtnp_(maxres6),vtnp_a(maxres6)
 314 c
 315       scale=.true.
 316       icount_scale=0
 317       if (lang.eq.1) then
 318         call sddir_precalc
 319       else if (lang.eq.2 .or. lang.eq.3) then
 320 #ifndef LANG0
 321         call stochastic_force(stochforcvec)
 322 #else
 323         write (iout,*)
 324      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 325 #ifdef MPI
 326         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 327 #endif
 328         stop
 329 #endif
 330       endif
 331       itime_scal=0
 332       do while (scale)
 333         icount_scale=icount_scale+1
 334         if (icount_scale.gt.maxcount_scale) then
 335           write (iout,*)
 336      &      "ERROR: too many attempts at scaling down the time step. ",
 337      &      "amax=",amax,"epdrift=",epdrift,
 338      &      "damax=",damax,"edriftmax=",edriftmax,
 339      &      "d_time=",d_time
 340           call flush(iout)
 341 #ifdef MPI
 342           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
 343 #endif
 344           stop
 345         endif
 346 c First step of the velocity Verlet algorithm
 347         if (lang.eq.2) then
 348 #ifndef LANG0
 349           call sd_verlet1
 350 #endif
 351         else if (lang.eq.3) then
 352 #ifndef LANG0
 353           call sd_verlet1_ciccotti
 354 #endif
 355         else if (lang.eq.1) then
 356           call sddir_verlet1
 357         else if (tnp1) then
 358           call tnp1_step1
 359         else if (tnp) then
 360           call tnp_step1
 361         else
 362           if (tnh) then
 363             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
 364             do i=0,2*nres
 365              do j=1,3
 366               d_t_old(j,i)=d_t_old(j,i)*scale_nh
 367              enddo
 368             enddo
 369           endif
 370           call verlet1
 371         endif
 372 c Build the chain from the newly calculated coordinates
 373         call chainbuild_cart
 374         if (rattle) call rattle1
 375         if (ntwe.ne.0) then
 376         if (large.and. mod(itime,ntwe).eq.0) then
 377           write (iout,*) "Cartesian and internal coordinates: step 1"
 378           call cartprint
 379           call intout
 380           write (iout,*) "dC"
 381           do i=0,nres
 382             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
 383      &      (dc(j,i+nres),j=1,3)
 384           enddo
 385           write (iout,*) "Accelerations"
 386           do i=0,nres
 387             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 388      &      (d_a(j,i+nres),j=1,3)
 389           enddo
 390           write (iout,*) "Velocities, step 1"
 391           do i=0,nres
 392             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 393      &      (d_t(j,i+nres),j=1,3)
 394           enddo
 395         endif
 396         endif
 397 #ifdef MPI
 398         tt0 = MPI_Wtime()
 399 #else
 400         tt0 = tcpu()
 401 #endif
 402 c Calculate energy and forces
 403         call zerograd
 404         call etotal(potEcomp)
 405 #ifdef TIMING_ENE
 406 #ifdef MPI
 407         t_etotal=t_etotal+MPI_Wtime()-tt0
 408 #else
 409         t_etotal=t_etotal+tcpu()-tt0
 410 #endif
 411 #endif
 412         E_old=potE
 413         potE=potEcomp(0)-potEcomp(20)
 414         call cartgrad
 415 c Get the new accelerations
 416         call lagrangian
 417 #ifdef MPI
 418         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 419 #else
 420         t_enegrad=t_enegrad+tcpu()-tt0
 421 #endif
 422 c Determine maximum acceleration and scale down the timestep if needed
 423         call max_accel
 424         amax=amax/(itime_scal+1)**2
 425         call predict_edrift(epdrift)
 426         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
 427 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 428           scale=.true.
 429           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 430      &      /dlog(2.0d0)+1
 431           itime_scal=itime_scal+ifac_time
 432 c          fac_time=dmin1(damax/amax,0.5d0)
 433           fac_time=0.5d0**ifac_time
 434           d_time=d_time*fac_time
 435           if (lang.eq.2 .or. lang.eq.3) then
 436 #ifndef LANG0
 437 c            write (iout,*) "Calling sd_verlet_setup: 1"
 438 c Rescale the stochastic forces and recalculate or restore
 439 c the matrices of tinker integrator
 440             if (itime_scal.gt.maxflag_stoch) then
 441               if (large) write (iout,'(a,i5,a)')
 442      &         "Calculate matrices for stochastic step;",
 443      &         " itime_scal ",itime_scal
 444               if (lang.eq.2) then
 445                 call sd_verlet_p_setup
 446               else
 447                 call sd_verlet_ciccotti_setup
 448               endif
 449               write (iout,'(2a,i3,a,i3,1h.)')
 450      &         "Warning: cannot store matrices for stochastic",
 451      &         " integration because the index",itime_scal,
 452      &         " is greater than",maxflag_stoch
 453               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 454      &         " integration Langevin algorithm for better efficiency."
 455             else if (flag_stoch(itime_scal)) then
 456               if (large) write (iout,'(a,i5,a,l1)')
 457      &         "Restore matrices for stochastic step; itime_scal ",
 458      &         itime_scal," flag ",flag_stoch(itime_scal)
 459               do i=1,dimen3
 460                 do j=1,dimen3
 461                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 462                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 463                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 464                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 465                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 466                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 467                 enddo
 468               enddo
 469             else
 470               if (large) write (iout,'(2a,i5,a,l1)')
 471      &         "Calculate & store matrices for stochastic step;",
 472      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 473               if (lang.eq.2) then
 474                 call sd_verlet_p_setup
 475               else
 476                 call sd_verlet_ciccotti_setup
 477               endif
 478               flag_stoch(ifac_time)=.true.
 479               do i=1,dimen3
 480                 do j=1,dimen3
 481                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 482                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 483                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 484                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 485                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 486                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 487                 enddo
 488               enddo
 489             endif
 490             fac_time=1.0d0/dsqrt(fac_time)
 491             do i=1,dimen3
 492               stochforcvec(i)=fac_time*stochforcvec(i)
 493             enddo
 494 #endif
 495           else if (lang.eq.1) then
 496 c Rescale the accelerations due to stochastic forces
 497             fac_time=1.0d0/dsqrt(fac_time)
 498             do i=1,dimen3
 499               d_as_work(i)=d_as_work(i)*fac_time
 500             enddo
 501           endif
 502           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 503      &      "itime",itime," Timestep scaled down to ",
 504      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 505         else
 506 c Second step of the velocity Verlet algorithm
 507           if (lang.eq.2) then
 508 #ifndef LANG0
 509             call sd_verlet2
 510 #endif
 511           else if (lang.eq.3) then
 512 #ifndef LANG0
 513             call sd_verlet2_ciccotti
 514 #endif
 515           else if (lang.eq.1) then
 516             call sddir_verlet2
 517           else if (tnp1) then
 518             call tnp1_step2
 519           else if (tnp) then
 520             call tnp_step2
 521           else
 522             call verlet2
 523             if (tnh) then
 524               call kinetic(EK)
 525               call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
 526               do i=0,2*nres
 527                do j=1,3
 528                 d_t(j,i)=d_t(j,i)*scale_nh
 529                enddo
 530               enddo
 531             endif
 532           endif
 533           if (rattle) call rattle2
 534           totT=totT+d_time
 535           if (d_time.ne.d_time0) then
 536             d_time=d_time0
 537 #ifndef   LANG0
 538             if (lang.eq.2 .or. lang.eq.3) then
 539               if (large) write (iout,'(a)')
 540      &         "Restore original matrices for stochastic step"
 541 c              write (iout,*) "Calling sd_verlet_setup: 2"
 542 c Restore the matrices of tinker integrator if the time step has been restored
 543               do i=1,dimen3
 544                 do j=1,dimen3
 545                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 546                   afric_mat(i,j)=afric0_mat(i,j,0)
 547                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 548                   prand_mat(i,j)=prand0_mat(i,j,0)
 549                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 550                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 551                 enddo
 552               enddo
 553             endif
 554 #endif
 555           endif
 556           scale=.false.
 557         endif
 558       enddo
 559 c Calculate the kinetic and the total energy and the kinetic temperature
 560       if (tnp .or. tnp1) then
 561        do i=0,2*nres
 562         do j=1,3
 563           d_t_old(j,i)=d_t(j,i)
 564           d_t(j,i)=d_t(j,i)/s_np
 565         enddo
 566        enddo
 567       endif
 568       call kinetic(EK)
 569       totE=EK+potE
 570 c diagnostics
 571 c      call kinetic1(EK1)
 572 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 573 c end diagnostics
 574 c Couple the system to Berendsen bath if needed
 575       if (tbf .and. lang.eq.0) then
 576         call verlet_bath
 577       endif
 578       kinetic_T=2.0d0/(dimen3*Rb)*EK
 579 c Backup the coordinates, velocities, and accelerations
 580       do i=0,2*nres
 581         do j=1,3
 582           dc_old(j,i)=dc(j,i)
 583           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 584           d_a_old(j,i)=d_a(j,i)
 585         enddo
 586       enddo
 587       if (ntwe.ne.0) then
 588       if (mod(itime,ntwe).eq.0) then
 589
 590        if(tnp .or. tnp1) then
 591         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
 592         H=(HNose1-H0)*s_np
 593 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
 594 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
 595 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 596           hhh=h
 597        endif
 598
 599        if(tnh) then
 600         HNose1=Hnose_nh(EK,potE)
 601         H=HNose1-H0
 602         hhh=h
 603 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 604        endif
 605
 606        if (large) then
 607         itnp=0
 608         do j=1,3
 609          itnp=itnp+1
 610          vtnp(itnp)=d_t(j,0)
 611         enddo
 612         do i=nnt,nct-1
 613          do j=1,3
 614           itnp=itnp+1
 615           vtnp(itnp)=d_t(j,i)
 616          enddo
 617         enddo
 618         do i=nnt,nct
 619          if (itype(i).ne.10) then
 620           inres=i+nres
 621           do j=1,3
 622            itnp=itnp+1
 623            vtnp(itnp)=d_t(j,inres)
 624           enddo
 625          endif
 626         enddo
 627
 628 c Transform velocities from UNRES coordinate space to cartesian and Gvec
 629 c eigenvector space
 630
 631         do i=1,dimen3
 632           vtnp_(i)=0.0d0
 633           vtnp_a(i)=0.0d0
 634           do j=1,dimen3
 635             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
 636             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
 637           enddo
 638           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
 639         enddo
 640
 641         do i=1,dimen3
 642          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
 643         enddo
 644
 645         write (iout,*) "Velocities, step 2"
 646         do i=0,nres
 647           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 648      &      (d_t(j,i+nres),j=1,3)
 649         enddo
 650        endif
 651       endif
 652       endif
 653       return
 654       end
 655 c-------------------------------------------------------------------------------
 656       subroutine RESPA_step(itime)
 657 c-------------------------------------------------------------------------------
 658 c  Perform a single RESPA step.
 659 c-------------------------------------------------------------------------------
 660       implicit real*8 (a-h,o-z)
 661       include 'DIMENSIONS'
 662 #ifdef MPI
 663       include 'mpif.h'
 664       integer IERROR,ERRCODE
 665 #endif
 666       include 'COMMON.SETUP'
 667       include 'COMMON.CONTROL'
 668       include 'COMMON.VAR'
 669       include 'COMMON.MD'
 670 #ifndef LANG0
 671       include 'COMMON.LANGEVIN'
 672 #else
 673       include 'COMMON.LANGEVIN.lang0'
 674 #endif
 675       include 'COMMON.CHAIN'
 676       include 'COMMON.DERIV'
 677       include 'COMMON.GEO'
 678       include 'COMMON.LOCAL'
 679       include 'COMMON.INTERACT'
 680       include 'COMMON.IOUNITS'
 681       include 'COMMON.NAMES'
 682       include 'COMMON.TIME1'
 683       double precision energia_short(0:n_ene),
 684      & energia_long(0:n_ene)
 685       double precision cm(3),L(3),vcm(3),incr(3)
 686       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
 687      & d_a_old0(3,0:maxres2)
 688       integer ilen,count,rstcount
 689       external ilen
 690       character*50 tytul
 691       integer maxcount_scale /10/
 692       common /gucio/ cm,energia_short
 693       double precision stochforcvec(MAXRES6)
 694       common /stochcalc/ stochforcvec
 695       integer itime
 696       logical scale
 697       double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
 698       common /cipiszcze/ itt
 699       itt=itime
 700       if (ntwe.ne.0) then
 701       if (large.and. mod(itime,ntwe).eq.0) then
 702         write (iout,*) "***************** RESPA itime",itime
 703         write (iout,*) "Cartesian and internal coordinates: step 0"
 704 c        call cartprint
 705         call pdbout(0.0d0,"cipiszcze",iout)
 706         call intout
 707         write (iout,*) "Accelerations from long-range forces"
 708         do i=0,nres
 709           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 710      &      (d_a(j,i+nres),j=1,3)
 711         enddo
 712         write (iout,*) "Velocities, step 0"
 713         do i=0,nres
 714           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 715      &      (d_t(j,i+nres),j=1,3)
 716         enddo
 717       endif
 718       endif
 719 c
 720 c Perform the initial RESPA step (increment velocities)
 721 c      write (iout,*) "*********************** RESPA ini"
 722       if (tnp1) then
 723           call tnp_respa_step1
 724       else if (tnp) then
 725           call tnp_respa_step1
 726       else
 727           if (tnh.and..not.xiresp) then
 728             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
 729             do i=0,2*nres
 730              do j=1,3
 731               d_t(j,i)=d_t(j,i)*scale_nh
 732              enddo
 733             enddo
 734           endif
 735           call RESPA_vel
 736       endif
 737
 738 cd       if(tnp .or. tnp1) then
 739 cd        write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
 740 cd       endif
 741
 742       if (ntwe.ne.0) then
 743       if (mod(itime,ntwe).eq.0 .and. large) then
 744         write (iout,*) "Velocities, end"
 745         do i=0,nres
 746           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 747      &      (d_t(j,i+nres),j=1,3)
 748         enddo
 749       endif
 750       endif
 751 c Compute the short-range forces
 752 #ifdef MPI
 753       tt0 =MPI_Wtime()
 754 #else
 755       tt0 = tcpu()
 756 #endif
 757 C 7/2/2009 commented out
 758 c      call zerograd
 759 c      call etotal_short(energia_short)
 760       if (tnp.or.tnp1) potE=energia_short(0)
 761 c      call cartgrad
 762 c      call lagrangian
 763 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
 764         do i=0,2*nres
 765           do j=1,3
 766             d_a(j,i)=d_a_short(j,i)
 767           enddo
 768         enddo
 769       if (ntwe.ne.0) then
 770       if (large.and. mod(itime,ntwe).eq.0) then
 771         write (iout,*) "energia_short",energia_short(0)
 772         write (iout,*) "Accelerations from short-range forces"
 773         do i=0,nres
 774           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 775      &      (d_a(j,i+nres),j=1,3)
 776         enddo
 777       endif
 778       endif
 779 #ifdef MPI
 780         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 781 #else
 782         t_enegrad=t_enegrad+tcpu()-tt0
 783 #endif
 784       do i=0,2*nres
 785         do j=1,3
 786           dc_old(j,i)=dc(j,i)
 787           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 788           d_a_old(j,i)=d_a(j,i)
 789         enddo
 790       enddo
 791 c 6/30/08 A-MTS: attempt at increasing the split number
 792       do i=0,2*nres
 793         do j=1,3
 794           dc_old0(j,i)=dc_old(j,i)
 795           d_t_old0(j,i)=d_t_old(j,i)
 796           d_a_old0(j,i)=d_a_old(j,i)
 797         enddo
 798       enddo
 799       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 800       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 801 c
 802       scale=.true.
 803       d_time0=d_time
 804       do while (scale)
 805
 806       scale=.false.
 807 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 808 c Split the time step
 809       d_time=d_time0/ntime_split
 810 c Perform the short-range RESPA steps (velocity Verlet increments of
 811 c positions and velocities using short-range forces)
 812 c      write (iout,*) "*********************** RESPA split"
 813       do itsplit=1,ntime_split
 814         if (lang.eq.1) then
 815           call sddir_precalc
 816         else if (lang.eq.2 .or. lang.eq.3) then
 817 #ifndef LANG0
 818           call stochastic_force(stochforcvec)
 819 #else
 820           write (iout,*)
 821      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 822 #ifdef MPI
 823           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 824 #endif
 825           stop
 826 #endif
 827         endif
 828 c First step of the velocity Verlet algorithm
 829         if (lang.eq.2) then
 830 #ifndef LANG0
 831           call sd_verlet1
 832 #endif
 833         else if (lang.eq.3) then
 834 #ifndef LANG0
 835           call sd_verlet1_ciccotti
 836 #endif
 837         else if (lang.eq.1) then
 838           call sddir_verlet1
 839         else if (tnp1) then
 840           call tnp1_respa_i_step1
 841         else if (tnp) then
 842           call tnp_respa_i_step1
 843         else
 844           if (tnh.and.xiresp) then
 845             call kinetic(EK)
 846             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
 847             do i=0,2*nres
 848              do j=1,3
 849               d_t_old(j,i)=d_t_old(j,i)*scale_nh
 850              enddo
 851             enddo
 852 cd            write(iout,*) "SSS1",itsplit,EK,scale_nh
 853           endif
 854           call verlet1
 855         endif
 856 c Build the chain from the newly calculated coordinates
 857         call chainbuild_cart
 858         if (rattle) call rattle1
 859         if (ntwe.ne.0) then
 860         if (large.and. mod(itime,ntwe).eq.0) then
 861           write (iout,*) "***** ITSPLIT",itsplit
 862           write (iout,*) "Cartesian and internal coordinates: step 1"
 863           call pdbout(0.0d0,"cipiszcze",iout)
 864 c          call cartprint
 865           call intout
 866           write (iout,*) "Velocities, step 1"
 867           do i=0,nres
 868             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 869      &        (d_t(j,i+nres),j=1,3)
 870           enddo
 871         endif
 872         endif
 873 #ifdef MPI
 874         tt0 = MPI_Wtime()
 875 #else
 876         tt0 = tcpu()
 877 #endif
 878 c Calculate energy and forces
 879         call zerograd
 880         call etotal_short(energia_short)
 881         E_old=potE
 882         potE=energia_short(0)
 883 #ifdef TIMING_ENE
 884 #ifdef MPI
 885         t_eshort=t_eshort+MPI_Wtime()-tt0
 886 #else
 887         t_eshort=t_eshort+tcpu()-tt0
 888 #endif
 889 #endif
 890         call cartgrad
 891 c Get the new accelerations
 892         call lagrangian
 893 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 894         do i=0,2*nres
 895           do j=1,3
 896             d_a_short(j,i)=d_a(j,i)
 897           enddo
 898         enddo
 899         if (ntwe.ne.0) then
 900         if (large.and. mod(itime,ntwe).eq.0) then
 901           write (iout,*)"energia_short",energia_short(0)
 902           write (iout,*) "Accelerations from short-range forces"
 903           do i=0,nres
 904             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 905      &        (d_a(j,i+nres),j=1,3)
 906           enddo
 907         endif
 908         endif
 909 c 6/30/08 A-MTS
 910 c Determine maximum acceleration and scale down the timestep if needed
 911         call max_accel
 912         amax=amax/ntime_split**2
 913         call predict_edrift(epdrift)
 914         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 915      &   write (iout,*) "amax",amax," damax",damax,
 916      &   " epdrift",epdrift," epdriftmax",epdriftmax
 917 c Exit loop and try with increased split number if the change of
 918 c acceleration is too big
 919         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 920           if (ntime_split.lt.maxtime_split) then
 921             scale=.true.
 922             ntime_split=ntime_split*2
 923             do i=0,2*nres
 924               do j=1,3
 925                 dc_old(j,i)=dc_old0(j,i)
 926                 d_t_old(j,i)=d_t_old0(j,i)
 927                 d_a_old(j,i)=d_a_old0(j,i)
 928               enddo
 929             enddo
 930             write (iout,*) "acceleration/energy drift too large",amax,
 931      &      epdrift," split increased to ",ntime_split," itime",itime,
 932      &       " itsplit",itsplit
 933             exit
 934           else
 935             write (iout,*)
 936      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 937      &       maxtime_split,
 938      &      " already reached!!! Trying to carry on!"
 939           endif
 940         endif
 941 #ifdef MPI
 942         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 943 #else
 944         t_enegrad=t_enegrad+tcpu()-tt0
 945 #endif
 946 c Second step of the velocity Verlet algorithm
 947         if (lang.eq.2) then
 948 #ifndef LANG0
 949           call sd_verlet2
 950 #endif
 951         else if (lang.eq.3) then
 952 #ifndef LANG0
 953           call sd_verlet2_ciccotti
 954 #endif
 955         else if (lang.eq.1) then
 956           call sddir_verlet2
 957         else if (tnp1) then
 958             call tnp1_respa_i_step2
 959         else if (tnp) then
 960             call tnp_respa_i_step2
 961         else
 962           call verlet2
 963           if (tnh.and.xiresp) then
 964             call kinetic(EK)
 965             call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
 966             do i=0,2*nres
 967              do j=1,3
 968               d_t(j,i)=d_t(j,i)*scale_nh
 969              enddo
 970             enddo
 971 cd            write(iout,*) "SSS2",itsplit,EK,scale_nh
 972           endif
 973         endif
 974         if (rattle) call rattle2
 975 c Backup the coordinates, velocities, and accelerations
 976         if (tnp .or. tnp1) then
 977          do i=0,2*nres
 978           do j=1,3
 979             d_t_old(j,i)=d_t(j,i)
 980             if (tnp) d_t(j,i)=d_t(j,i)/s_np
 981             if (tnp1) d_t(j,i)=d_t(j,i)/s_np
 982           enddo
 983          enddo
 984         endif
 985
 986         do i=0,2*nres
 987           do j=1,3
 988             dc_old(j,i)=dc(j,i)
 989             if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 990             d_a_old(j,i)=d_a(j,i)
 991           enddo
 992         enddo
 993       enddo
 994
 995       enddo ! while scale
 996
 997 c Restore the time step
 998       d_time=d_time0
 999 c Compute long-range forces
1000 #ifdef MPI
1001       tt0 =MPI_Wtime()
1002 #else
1003       tt0 = tcpu()
1004 #endif
1005       call zerograd
1006       call etotal_long(energia_long)
1007       E_long=energia_long(0)
1008       potE=energia_short(0)+energia_long(0)
1009 #ifdef TIMING_ENE
1010 #ifdef MPI
1011         t_elong=t_elong+MPI_Wtime()-tt0
1012 #else
1013         t_elong=t_elong+tcpu()-tt0
1014 #endif
1015 #endif
1016       call cartgrad
1017       call lagrangian
1018 #ifdef MPI
1019         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1020 #else
1021         t_enegrad=t_enegrad+tcpu()-tt0
1022 #endif
1023 c Compute accelerations from long-range forces
1024       if (ntwe.ne.0) then
1025       if (large.and. mod(itime,ntwe).eq.0) then
1026         write (iout,*) "energia_long",energia_long(0)
1027         write (iout,*) "Cartesian and internal coordinates: step 2"
1028 c        call cartprint
1029         call pdbout(0.0d0,"cipiszcze",iout)
1030         call intout
1031         write (iout,*) "Accelerations from long-range forces"
1032         do i=0,nres
1033           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1034      &      (d_a(j,i+nres),j=1,3)
1035         enddo
1036         write (iout,*) "Velocities, step 2"
1037         do i=0,nres
1038           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1039      &      (d_t(j,i+nres),j=1,3)
1040         enddo
1041       endif
1042       endif
1043 c Compute the final RESPA step (increment velocities)
1044 c      write (iout,*) "*********************** RESPA fin"
1045       if (tnp1) then
1046           call tnp_respa_step2
1047       else if (tnp) then
1048           call tnp_respa_step2
1049       else
1050           call RESPA_vel
1051           if (tnh.and..not.xiresp) then
1052             call kinetic(EK)
1053             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
1054             do i=0,2*nres
1055              do j=1,3
1056               d_t(j,i)=d_t(j,i)*scale_nh
1057              enddo
1058             enddo
1059           endif
1060       endif
1061
1062         if (tnp .or. tnp1) then
1063          do i=0,2*nres
1064           do j=1,3
1065             d_t(j,i)=d_t_old(j,i)/s_np
1066           enddo
1067          enddo
1068         endif
1069
1070 c Compute the complete potential energy
1071       do i=0,n_ene
1072         potEcomp(i)=energia_short(i)+energia_long(i)
1073       enddo
1074       potE=potEcomp(0)-potEcomp(20)
1075 c      potE=energia_short(0)+energia_long(0)
1076       totT=totT+d_time
1077 c Calculate the kinetic and the total energy and the kinetic temperature
1078       call kinetic(EK)
1079       totE=EK+potE
1080 c Couple the system to Berendsen bath if needed
1081       if (tbf .and. lang.eq.0) then
1082         call verlet_bath
1083       endif
1084       kinetic_T=2.0d0/(dimen3*Rb)*EK
1085 c Backup the coordinates, velocities, and accelerations
1086       if (ntwe.ne.0) then
1087       if (mod(itime,ntwe).eq.0 .and. large) then
1088         write (iout,*) "Velocities, end"
1089         do i=0,nres
1090           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1091      &      (d_t(j,i+nres),j=1,3)
1092         enddo
1093       endif
1094
1095       if (mod(itime,ntwe).eq.0) then
1096
1097        if(tnp .or. tnp1) then
1098 #ifndef G77
1099         write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
1100      &                          E_long,energia_short(0)
1101 #else
1102         write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
1103      &                          E_long,energia_short(0)
1104 #endif
1105         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1106         H=(HNose1-H0)*s_np
1107 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
1108 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
1109 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1110           hhh=h
1111 cd        write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
1112        endif
1113
1114        if(tnh) then
1115         HNose1=Hnose_nh(EK,potE)
1116         H=HNose1-H0
1117 cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
1118         hhh=h
1119        endif
1120
1121
1122        if (large) then
1123        itnp=0
1124        do j=1,3
1125         itnp=itnp+1
1126         vtnp(itnp)=d_t(j,0)
1127        enddo
1128        do i=nnt,nct-1
1129         do j=1,3
1130           itnp=itnp+1
1131           vtnp(itnp)=d_t(j,i)
1132         enddo
1133        enddo
1134        do i=nnt,nct
1135         if (itype(i).ne.10) then
1136           inres=i+nres
1137           do j=1,3
1138            itnp=itnp+1
1139            vtnp(itnp)=d_t(j,inres)
1140           enddo
1141         endif
1142        enddo
1143
1144 c Transform velocities from UNRES coordinate space to cartesian and Gvec
1145 c eigenvector space
1146
1147         do i=1,dimen3
1148           vtnp_(i)=0.0d0
1149           vtnp_a(i)=0.0d0
1150           do j=1,dimen3
1151             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
1152             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
1153           enddo
1154           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
1155         enddo
1156
1157         do i=1,dimen3
1158          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
1159         enddo
1160
1161        endif
1162       endif
1163       endif
1164
1165
1166       return
1167       end
1168 c---------------------------------------------------------------------
1169       subroutine RESPA_vel
1170 c  First and last RESPA step (incrementing velocities using long-range
1171 c  forces).
1172       implicit real*8 (a-h,o-z)
1173       include 'DIMENSIONS'
1174       include 'COMMON.CONTROL'
1175       include 'COMMON.VAR'
1176       include 'COMMON.MD'
1177       include 'COMMON.CHAIN'
1178       include 'COMMON.DERIV'
1179       include 'COMMON.GEO'
1180       include 'COMMON.LOCAL'
1181       include 'COMMON.INTERACT'
1182       include 'COMMON.IOUNITS'
1183       include 'COMMON.NAMES'
1184       do j=1,3
1185         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1186       enddo
1187       do i=nnt,nct-1
1188         do j=1,3
1189           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1190         enddo
1191       enddo
1192       do i=nnt,nct
1193         if (itype(i).ne.10) then
1194           inres=i+nres
1195           do j=1,3
1196             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1197           enddo
1198         endif
1199       enddo
1200       return
1201       end
1202 c-----------------------------------------------------------------
1203       subroutine verlet1
1204 c Applying velocity Verlet algorithm - step 1 to coordinates
1205       implicit real*8 (a-h,o-z)
1206       include 'DIMENSIONS'
1207       include 'COMMON.CONTROL'
1208       include 'COMMON.VAR'
1209       include 'COMMON.MD'
1210       include 'COMMON.CHAIN'
1211       include 'COMMON.DERIV'
1212       include 'COMMON.GEO'
1213       include 'COMMON.LOCAL'
1214       include 'COMMON.INTERACT'
1215       include 'COMMON.IOUNITS'
1216       include 'COMMON.NAMES'
1217       double precision adt,adt2
1218
1219 #ifdef DEBUG
1220       write (iout,*) "VELVERLET1 START: DC"
1221       do i=0,nres
1222         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1223      &   (dc(j,i+nres),j=1,3)
1224       enddo
1225 #endif
1226       do j=1,3
1227         adt=d_a_old(j,0)*d_time
1228         adt2=0.5d0*adt
1229         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1230         d_t_new(j,0)=d_t_old(j,0)+adt2
1231         d_t(j,0)=d_t_old(j,0)+adt
1232       enddo
1233       do i=nnt,nct-1
1234         do j=1,3
1235           adt=d_a_old(j,i)*d_time
1236           adt2=0.5d0*adt
1237           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1238           d_t_new(j,i)=d_t_old(j,i)+adt2
1239           d_t(j,i)=d_t_old(j,i)+adt
1240         enddo
1241       enddo
1242       do i=nnt,nct
1243         if (itype(i).ne.10) then
1244           inres=i+nres
1245           do j=1,3
1246             adt=d_a_old(j,inres)*d_time
1247             adt2=0.5d0*adt
1248             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1249             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1250             d_t(j,inres)=d_t_old(j,inres)+adt
1251           enddo
1252         endif
1253       enddo
1254 #ifdef DEBUG
1255       write (iout,*) "VELVERLET1 END: DC"
1256       do i=0,nres
1257         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1258      &   (dc(j,i+nres),j=1,3)
1259       enddo
1260 #endif
1261       return
1262       end
1263 c---------------------------------------------------------------------
1264       subroutine verlet2
1265 c  Step 2 of the velocity Verlet algorithm: update velocities
1266       implicit real*8 (a-h,o-z)
1267       include 'DIMENSIONS'
1268       include 'COMMON.CONTROL'
1269       include 'COMMON.VAR'
1270       include 'COMMON.MD'
1271       include 'COMMON.CHAIN'
1272       include 'COMMON.DERIV'
1273       include 'COMMON.GEO'
1274       include 'COMMON.LOCAL'
1275       include 'COMMON.INTERACT'
1276       include 'COMMON.IOUNITS'
1277       include 'COMMON.NAMES'
1278       do j=1,3
1279         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1280       enddo
1281       do i=nnt,nct-1
1282         do j=1,3
1283           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1284         enddo
1285       enddo
1286       do i=nnt,nct
1287         if (itype(i).ne.10) then
1288           inres=i+nres
1289           do j=1,3
1290             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1291           enddo
1292         endif
1293       enddo
1294       return
1295       end
1296 c-----------------------------------------------------------------
1297       subroutine sddir_precalc
1298 c Applying velocity Verlet algorithm - step 1 to coordinates
1299       implicit real*8 (a-h,o-z)
1300       include 'DIMENSIONS'
1301 #ifdef MPI
1302       include 'mpif.h'
1303 #endif
1304       include 'COMMON.CONTROL'
1305       include 'COMMON.VAR'
1306       include 'COMMON.MD'
1307 #ifndef LANG0
1308       include 'COMMON.LANGEVIN'
1309 #else
1310       include 'COMMON.LANGEVIN.lang0'
1311 #endif
1312       include 'COMMON.CHAIN'
1313       include 'COMMON.DERIV'
1314       include 'COMMON.GEO'
1315       include 'COMMON.LOCAL'
1316       include 'COMMON.INTERACT'
1317       include 'COMMON.IOUNITS'
1318       include 'COMMON.NAMES'
1319       include 'COMMON.TIME1'
1320       double precision stochforcvec(MAXRES6)
1321       common /stochcalc/ stochforcvec
1322 c
1323 c Compute friction and stochastic forces
1324 c
1325 #ifdef MPI
1326       time00=MPI_Wtime()
1327 #else
1328       time00=tcpu()
1329 #endif
1330       call friction_force
1331 #ifdef MPI
1332       time_fric=time_fric+MPI_Wtime()-time00
1333       time00=MPI_Wtime()
1334 #else
1335       time_fric=time_fric+tcpu()-time00
1336       time00=tcpu()
1337 #endif
1338       call stochastic_force(stochforcvec)
1339 #ifdef MPI
1340       time_stoch=time_stoch+MPI_Wtime()-time00
1341 #else
1342       time_stoch=time_stoch+tcpu()-time00
1343 #endif
1344 c
1345 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1346 c forces (d_as_work)
1347 c
1348       call ginv_mult(fric_work, d_af_work)
1349       call ginv_mult(stochforcvec, d_as_work)
1350       return
1351       end
1352 c---------------------------------------------------------------------
1353       subroutine sddir_verlet1
1354 c Applying velocity Verlet algorithm - step 1 to velocities
1355       implicit real*8 (a-h,o-z)
1356       include 'DIMENSIONS'
1357       include 'COMMON.CONTROL'
1358       include 'COMMON.VAR'
1359       include 'COMMON.MD'
1360 #ifndef LANG0
1361       include 'COMMON.LANGEVIN'
1362 #else
1363       include 'COMMON.LANGEVIN.lang0'
1364 #endif
1365       include 'COMMON.CHAIN'
1366       include 'COMMON.DERIV'
1367       include 'COMMON.GEO'
1368       include 'COMMON.LOCAL'
1369       include 'COMMON.INTERACT'
1370       include 'COMMON.IOUNITS'
1371       include 'COMMON.NAMES'
1372 c Revised 3/31/05 AL: correlation between random contributions to
1373 c position and velocity increments included.
1374       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1375       double precision adt,adt2
1376 c
1377 c Add the contribution from BOTH friction and stochastic force to the
1378 c coordinates, but ONLY the contribution from the friction forces to velocities
1379 c
1380       do j=1,3
1381         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1382         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1383         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1384         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1385         d_t(j,0)=d_t_old(j,0)+adt
1386       enddo
1387       ind=3
1388       do i=nnt,nct-1
1389         do j=1,3
1390           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1391           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1392           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1393           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1394           d_t(j,i)=d_t_old(j,i)+adt
1395         enddo
1396         ind=ind+3
1397       enddo
1398       do i=nnt,nct
1399         if (itype(i).ne.10) then
1400           inres=i+nres
1401           do j=1,3
1402             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1403             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1404             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1405             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1406             d_t(j,inres)=d_t_old(j,inres)+adt
1407           enddo
1408           ind=ind+3
1409         endif
1410       enddo
1411       return
1412       end
1413 c---------------------------------------------------------------------
1414       subroutine sddir_verlet2
1415 c  Calculating the adjusted velocities for accelerations
1416       implicit real*8 (a-h,o-z)
1417       include 'DIMENSIONS'
1418       include 'COMMON.CONTROL'
1419       include 'COMMON.VAR'
1420       include 'COMMON.MD'
1421 #ifndef LANG0
1422       include 'COMMON.LANGEVIN'
1423 #else
1424       include 'COMMON.LANGEVIN.lang0'
1425 #endif
1426       include 'COMMON.CHAIN'
1427       include 'COMMON.DERIV'
1428       include 'COMMON.GEO'
1429       include 'COMMON.LOCAL'
1430       include 'COMMON.INTERACT'
1431       include 'COMMON.IOUNITS'
1432       include 'COMMON.NAMES'
1433       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1434       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1435 c Revised 3/31/05 AL: correlation between random contributions to
1436 c position and velocity increments included.
1437 c The correlation coefficients are calculated at low-friction limit.
1438 c Also, friction forces are now not calculated with new velocities.
1439
1440 c      call friction_force
1441       call stochastic_force(stochforcvec)
1442 c
1443 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1444 c forces (d_as_work)
1445 c
1446       call ginv_mult(stochforcvec, d_as_work1)
1447
1448 c
1449 c Update velocities
1450 c
1451       do j=1,3
1452         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1453      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1454       enddo
1455       ind=3
1456       do i=nnt,nct-1
1457         do j=1,3
1458           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1459      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1460         enddo
1461         ind=ind+3
1462       enddo
1463       do i=nnt,nct
1464         if (itype(i).ne.10) then
1465           inres=i+nres
1466           do j=1,3
1467             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1468      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1469      &       +cos60*d_as_work1(ind+j))*d_time
1470           enddo
1471           ind=ind+3
1472         endif
1473       enddo
1474       return
1475       end
1476 c---------------------------------------------------------------------
1477       subroutine max_accel
1478 c
1479 c Find the maximum difference in the accelerations of the the sites
1480 c at the beginning and the end of the time step.
1481 c
1482       implicit real*8 (a-h,o-z)
1483       include 'DIMENSIONS'
1484       include 'COMMON.CONTROL'
1485       include 'COMMON.VAR'
1486       include 'COMMON.MD'
1487       include 'COMMON.CHAIN'
1488       include 'COMMON.DERIV'
1489       include 'COMMON.GEO'
1490       include 'COMMON.LOCAL'
1491       include 'COMMON.INTERACT'
1492       include 'COMMON.IOUNITS'
1493       double precision aux(3),accel(3),accel_old(3),dacc
1494       do j=1,3
1495 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1496          accel_old(j)=d_a_old(j,0)
1497          accel(j)=d_a(j,0)
1498       enddo
1499       amax=0.0d0
1500       do i=nnt,nct
1501 c Backbone
1502         if (i.lt.nct) then
1503 c 7/3/08 changed to asymmetric difference
1504           do j=1,3
1505 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1506             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1507             accel(j)=accel(j)+0.5d0*d_a(j,i)
1508 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1509             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1510               dacc=dabs(accel(j)-accel_old(j))
1511               if (dacc.gt.amax) amax=dacc
1512             endif
1513           enddo
1514         endif
1515       enddo
1516 c Side chains
1517       do j=1,3
1518 c        accel(j)=aux(j)
1519         accel_old(j)=d_a_old(j,0)
1520         accel(j)=d_a(j,0)
1521       enddo
1522       if (nnt.eq.2) then
1523         do j=1,3
1524           accel_old(j)=accel_old(j)+d_a_old(j,1)
1525           accel(j)=accel(j)+d_a(j,1)
1526         enddo
1527       endif
1528       do i=nnt,nct
1529         if (itype(i).ne.10) then
1530           do j=1,3
1531 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1532             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1533             accel(j)=accel(j)+d_a(j,i+nres)
1534           enddo
1535         endif
1536         do j=1,3
1537 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1538           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1539             dacc=dabs(accel(j)-accel_old(j))
1540             if (dacc.gt.amax) amax=dacc
1541           endif
1542         enddo
1543         do j=1,3
1544           accel_old(j)=accel_old(j)+d_a_old(j,i)
1545           accel(j)=accel(j)+d_a(j,i)
1546 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1547         enddo
1548       enddo
1549       return
1550       end
1551 c---------------------------------------------------------------------
1552       subroutine predict_edrift(epdrift)
1553 c
1554 c Predict the drift of the potential energy
1555 c
1556       implicit real*8 (a-h,o-z)
1557       include 'DIMENSIONS'
1558       include 'COMMON.CONTROL'
1559       include 'COMMON.VAR'
1560       include 'COMMON.MD'
1561       include 'COMMON.CHAIN'
1562       include 'COMMON.DERIV'
1563       include 'COMMON.GEO'
1564       include 'COMMON.LOCAL'
1565       include 'COMMON.INTERACT'
1566       include 'COMMON.IOUNITS'
1567       include 'COMMON.MUCA'
1568       double precision epdrift,epdriftij
1569 c Drift of the potential energy
1570       epdrift=0.0d0
1571       do i=nnt,nct
1572 c Backbone
1573         if (i.lt.nct) then
1574           do j=1,3
1575             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1576             if (lmuca) epdriftij=epdriftij*factor
1577 c            write (iout,*) "back",i,j,epdriftij
1578             if (epdriftij.gt.epdrift) epdrift=epdriftij
1579           enddo
1580         endif
1581 c Side chains
1582         if (itype(i).ne.10) then
1583           do j=1,3
1584             epdriftij=
1585      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1586             if (lmuca) epdriftij=epdriftij*factor
1587 c            write (iout,*) "side",i,j,epdriftij
1588             if (epdriftij.gt.epdrift) epdrift=epdriftij
1589           enddo
1590         endif
1591       enddo
1592       epdrift=0.5d0*epdrift*d_time*d_time
1593 c      write (iout,*) "epdrift",epdrift
1594       return
1595       end
1596 c-----------------------------------------------------------------------
1597       subroutine verlet_bath
1598 c
1599 c  Coupling to the thermostat by using the Berendsen algorithm
1600 c
1601       implicit real*8 (a-h,o-z)
1602       include 'DIMENSIONS'
1603       include 'COMMON.CONTROL'
1604       include 'COMMON.VAR'
1605       include 'COMMON.MD'
1606       include 'COMMON.CHAIN'
1607       include 'COMMON.DERIV'
1608       include 'COMMON.GEO'
1609       include 'COMMON.LOCAL'
1610       include 'COMMON.INTERACT'
1611       include 'COMMON.IOUNITS'
1612       include 'COMMON.NAMES'
1613       double precision T_half,fact
1614 c
1615       T_half=2.0d0/(dimen3*Rb)*EK
1616       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1617 c      write(iout,*) "T_half", T_half
1618 c      write(iout,*) "EK", EK
1619 c      write(iout,*) "fact", fact
1620       do j=1,3
1621         d_t(j,0)=fact*d_t(j,0)
1622       enddo
1623       do i=nnt,nct-1
1624         do j=1,3
1625           d_t(j,i)=fact*d_t(j,i)
1626         enddo
1627       enddo
1628       do i=nnt,nct
1629         if (itype(i).ne.10) then
1630           inres=i+nres
1631           do j=1,3
1632             d_t(j,inres)=fact*d_t(j,inres)
1633           enddo
1634         endif
1635       enddo
1636       return
1637       end
1638 c---------------------------------------------------------
1639       subroutine init_MD
1640 c  Set up the initial conditions of a MD simulation
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643 #ifdef MP
1644       include 'mpif.h'
1645       character*16 form
1646       integer IERROR,ERRCODE
1647 #endif
1648       include 'COMMON.SETUP'
1649       include 'COMMON.CONTROL'
1650       include 'COMMON.VAR'
1651       include 'COMMON.MD'
1652 #ifndef LANG0
1653       include 'COMMON.LANGEVIN'
1654 #else
1655       include 'COMMON.LANGEVIN.lang0'
1656 #endif
1657       include 'COMMON.CHAIN'
1658       include 'COMMON.DERIV'
1659       include 'COMMON.GEO'
1660       include 'COMMON.LOCAL'
1661       include 'COMMON.INTERACT'
1662       include 'COMMON.IOUNITS'
1663       include 'COMMON.NAMES'
1664       include 'COMMON.REMD'
1665       real*8 energia_long(0:n_ene),
1666      &  energia_short(0:n_ene),vcm(3),incr(3),E_short
1667       double precision cm(3),L(3),xv,sigv,lowb,highb
1668       double precision varia(maxvar)
1669       character*256 qstr
1670       integer ilen
1671       external ilen
1672       character*50 tytul
1673       logical file_exist
1674       common /gucio/ cm
1675       d_time0=d_time
1676 c      write(iout,*) "d_time", d_time
1677 c Compute the standard deviations of stochastic forces for Langevin dynamics
1678 c if the friction coefficients do not depend on surface area
1679       if (lang.gt.0 .and. .not.surfarea) then
1680         do i=nnt,nct-1
1681           stdforcp(i)=stdfp*dsqrt(gamp)
1682         enddo
1683         do i=nnt,nct
1684           stdforcsc(i)=stdfsc(iabs(itype(i)))
1685      &                *dsqrt(gamsc(iabs(itype(i))))
1686         enddo
1687       endif
1688 c Open the pdb file for snapshotshots
1689 #ifdef MPI
1690       if(mdpdb) then
1691         if (ilen(tmpdir).gt.0)
1692      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1693      &      liczba(:ilen(liczba))//".pdb")
1694         open(ipdb,
1695      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1696      &  //".pdb")
1697       else
1698 #ifdef NOXDR
1699         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1700      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1701      &      liczba(:ilen(liczba))//".x")
1702         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1703      &  //".x"
1704 #else
1705         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1706      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1707      &      liczba(:ilen(liczba))//".cx")
1708         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1709      &  //".cx"
1710 #endif
1711       endif
1712 #else
1713       if(mdpdb) then
1714          if (ilen(tmpdir).gt.0)
1715      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1716          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1717       else
1718          if (ilen(tmpdir).gt.0)
1719      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1720          cartname=prefix(:ilen(prefix))//"_MD.cx"
1721       endif
1722 #endif
1723       if (usampl) then
1724         write (qstr,'(256(1h ))')
1725         ipos=1
1726         do i=1,nfrag
1727           iq = qinfrag(i,iset)*10
1728           iw = wfrag(i,iset)/100
1729           if (iw.gt.0) then
1730             if(me.eq.king.or..not.out1file)
1731      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1732             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1733             ipos=ipos+7
1734           endif
1735         enddo
1736         do i=1,npair
1737           iq = qinpair(i,iset)*10
1738           iw = wpair(i,iset)/100
1739           if (iw.gt.0) then
1740             if(me.eq.king.or..not.out1file)
1741      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1742             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1743             ipos=ipos+7
1744           endif
1745         enddo
1746 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1747 #ifdef NOXDR
1748 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1749 #else
1750 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1751 #endif
1752 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1753       endif
1754       icg=1
1755       if (rest) then
1756        if (restart1file) then
1757          if (me.eq.king)
1758      &     inquire(file=mremd_rst_name,exist=file_exist)
1759            write (*,*) me," Before broadcast: file_exist",file_exist
1760 #ifdef MPI
1761          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1762      &          IERR)
1763          write (*,*) me," After broadcast: file_exist",file_exist
1764 #endif
1765 c        inquire(file=mremd_rst_name,exist=file_exist)
1766         if(me.eq.king.or..not.out1file)
1767      &   write(iout,*) "Initial state read by master and distributed"
1768        else
1769          if (ilen(tmpdir).gt.0)
1770      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1771      &      //liczba(:ilen(liczba))//'.rst')
1772         inquire(file=rest2name,exist=file_exist)
1773        endif
1774        if(file_exist) then
1775          if(.not.restart1file) then
1776            if(me.eq.king.or..not.out1file)
1777      &      write(iout,*) "Initial state will be read from file ",
1778      &      rest2name(:ilen(rest2name))
1779            call readrst
1780          endif
1781          call rescale_weights(t_bath)
1782        else
1783         if(me.eq.king.or..not.out1file)then
1784          if (restart1file) then
1785           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1786      &       " does not exist"
1787          else
1788           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1789      &       " does not exist"
1790          endif
1791          write(iout,*) "Initial velocities randomly generated"
1792         endif
1793         call random_vel
1794         totT=0.0d0
1795        endif
1796       else
1797 c Generate initial velocities
1798         if(me.eq.king.or..not.out1file)
1799      &   write(iout,*) "Initial velocities randomly generated"
1800         call random_vel
1801         totT=0.0d0
1802       endif
1803 c      rest2name = prefix(:ilen(prefix))//'.rst'
1804       if(me.eq.king.or..not.out1file)then
1805        write (iout,*) "Initial velocities"
1806        do i=0,nres
1807          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1808      &   (d_t(j,i+nres),j=1,3)
1809        enddo
1810        call flush(iout)
1811 c  Zeroing the total angular momentum of the system
1812        write(iout,*) "Calling the zero-angular
1813      & momentum subroutine"
1814       endif
1815       call inertia_tensor
1816 c  Getting the potential energy and forces and velocities and accelerations
1817       call vcm_vel(vcm)
1818 c      write (iout,*) "velocity of the center of the mass:"
1819 c      write (iout,*) (vcm(j),j=1,3)
1820       do j=1,3
1821         d_t(j,0)=d_t(j,0)-vcm(j)
1822       enddo
1823 c Removing the velocity of the center of mass
1824       call vcm_vel(vcm)
1825       if(me.eq.king.or..not.out1file)then
1826        write (iout,*) "vcm right after adjustment:"
1827        write (iout,*) (vcm(j),j=1,3)
1828        call flush(iout)
1829       endif
1830       if (.not.rest) then
1831          call chainbuild
1832          if(iranconf.ne.0) then
1833           if (overlapsc) then
1834            print *, 'Calling OVERLAP_SC'
1835            call overlap_sc(fail)
1836           endif
1837
1838           if (searchsc) then
1839            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1840            print *,'SC_move',nft_sc,etot
1841            if(me.eq.king.or..not.out1file)
1842      &      write(iout,*) 'SC_move',nft_sc,etot
1843           endif
1844
1845           if(dccart)then
1846            print *, 'Calling MINIM_DC'
1847            call minim_dc(etot,iretcode,nfun)
1848           else
1849            call geom_to_var(nvar,varia)
1850            print *,'Calling MINIMIZE.'
1851            call minimize(etot,varia,iretcode,nfun)
1852            call var_to_geom(nvar,varia)
1853           endif
1854           if(me.eq.king.or..not.out1file)
1855      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1856          endif
1857       endif
1858       call chainbuild_cart
1859       call kinetic(EK)
1860       if (tbf) then
1861         call verlet_bath(EK)
1862       endif
1863       kinetic_T=2.0d0/(dimen3*Rb)*EK
1864       if(me.eq.king.or..not.out1file)then
1865        call cartprint
1866        call intout
1867       endif
1868 #ifdef MPI
1869       tt0=MPI_Wtime()
1870 #else
1871       tt0=tcpu()
1872 #endif
1873       call zerograd
1874       call etotal(potEcomp)
1875 #ifdef TIMING_ENE
1876 #ifdef MPI
1877       t_etotal=t_etotal+MPI_Wtime()-tt0
1878 #else
1879       t_etotal=t_etotal+tcpu()-tt0
1880 #endif
1881 #endif
1882       potE=potEcomp(0)
1883
1884       if(tnp .or. tnp1) then
1885        s_np=1.0
1886        pi_np=0.0
1887        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1888        H0=Hnose1
1889        write(iout,*) 'H0= ',H0
1890       endif
1891
1892       if(tnh) then
1893        HNose1=Hnose_nh(EK,potE)
1894        H0=HNose1
1895        write (iout,*) 'H0= ',H0
1896       endif
1897
1898       if (hmc.gt.0) then
1899          hmc_acc=0
1900          hmc_etot=potE+EK
1901           if(me.eq.king.or..not.out1file)
1902      &       write(iout,*) 'HMC',hmc_etot,potE,EK
1903          do i=1,2*nres
1904            do j=1,3
1905             dc_hmc(j,i)=dc(j,i)
1906            enddo
1907          enddo
1908       endif
1909
1910       call cartgrad
1911       call lagrangian
1912       call max_accel
1913       if (amax*d_time .gt. dvmax) then
1914         d_time=d_time*dvmax/amax
1915         if(me.eq.king.or..not.out1file) write (iout,*)
1916      &   "Time step reduced to",d_time,
1917      &   " because of too large initial acceleration."
1918       endif
1919       if(me.eq.king.or..not.out1file)then
1920        write(iout,*) "Potential energy and its components"
1921        call enerprint(potEcomp)
1922 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1923       endif
1924       potE=potEcomp(0)-potEcomp(20)
1925       totE=EK+potE
1926       itime=0
1927       if (ntwe.ne.0) call statout(itime)
1928       if(me.eq.king.or..not.out1file)
1929      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1930      &   " Kinetic energy",EK," potential energy",potE,
1931      &   " total energy",totE," maximum acceleration ",
1932      &   amax
1933       if (large) then
1934         write (iout,*) "Initial coordinates"
1935         do i=1,nres
1936           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1937      &    (c(j,i+nres),j=1,3)
1938         enddo
1939         write (iout,*) "Initial dC"
1940         do i=0,nres
1941           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1942      &    (dc(j,i+nres),j=1,3)
1943         enddo
1944         write (iout,*) "Initial velocities"
1945         write (iout,"(13x,' backbone ',23x,' side chain')")
1946         do i=0,nres
1947           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1948      &    (d_t(j,i+nres),j=1,3)
1949         enddo
1950         write (iout,*) "Initial accelerations"
1951         do i=0,nres
1952 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1953           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1954      &    (d_a(j,i+nres),j=1,3)
1955         enddo
1956       endif
1957       do i=0,2*nres
1958         do j=1,3
1959           dc_old(j,i)=dc(j,i)
1960           d_t_old(j,i)=d_t(j,i)
1961           d_a_old(j,i)=d_a(j,i)
1962         enddo
1963 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1964       enddo
1965       if (RESPA) then
1966 #ifdef MPI
1967       tt0 =MPI_Wtime()
1968 #else
1969       tt0 = tcpu()
1970 #endif
1971         call zerograd
1972         call etotal_short(energia_short)
1973 #ifdef TIMING_ENE
1974 #ifdef MPI
1975         t_eshort=t_eshort+MPI_Wtime()-tt0
1976 #else
1977         t_eshort=t_eshort+tcpu()-tt0
1978 #endif
1979 #endif
1980
1981         if(tnp .or. tnp1) then
1982          E_short=energia_short(0)
1983          HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1984          Csplit=Hnose1
1985 c         Csplit =110
1986 c_new_var_csplit          Csplit=H0-E_long
1987 c          Csplit = H0-energia_short(0)
1988           write(iout,*) 'Csplit= ',Csplit
1989         endif
1990
1991
1992         call cartgrad
1993         call lagrangian
1994         if(.not.out1file .and. large) then
1995           write (iout,*) "energia_long",energia_long(0),
1996      &     " energia_short",energia_short(0),
1997      &     " total",energia_long(0)+energia_short(0)
1998           write (iout,*) "Initial fast-force accelerations"
1999           do i=0,nres
2000             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2001      &      (d_a(j,i+nres),j=1,3)
2002           enddo
2003         endif
2004 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2005         do i=0,2*nres
2006           do j=1,3
2007             d_a_short(j,i)=d_a(j,i)
2008           enddo
2009         enddo
2010 #ifdef MPI
2011         tt0=MPI_Wtime()
2012 #else
2013         tt0=tcpu()
2014 #endif
2015         call zerograd
2016         call etotal_long(energia_long)
2017 #ifdef TIMING_ENE
2018 #ifdef MPI
2019         t_elong=t_elong+MPI_Wtime()-tt0
2020 #else
2021         t_elong=t_elong+tcpu()-tt0
2022 #endif
2023 #endif
2024         call cartgrad
2025         call lagrangian
2026         if(.not.out1file .and. large) then
2027           write (iout,*) "energia_long",energia_long(0)
2028           write (iout,*) "Initial slow-force accelerations"
2029           do i=0,nres
2030             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2031      &      (d_a(j,i+nres),j=1,3)
2032           enddo
2033         endif
2034 #ifdef MPI
2035         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2036 #else
2037         t_enegrad=t_enegrad+tcpu()-tt0
2038 #endif
2039       endif
2040
2041
2042
2043       return
2044       end
2045 c-----------------------------------------------------------
2046       subroutine random_vel
2047       implicit real*8 (a-h,o-z)
2048       include 'DIMENSIONS'
2049       include 'COMMON.CONTROL'
2050       include 'COMMON.VAR'
2051       include 'COMMON.MD'
2052 #ifndef LANG0
2053       include 'COMMON.LANGEVIN'
2054 #else
2055       include 'COMMON.LANGEVIN.lang0'
2056 #endif
2057       include 'COMMON.CHAIN'
2058       include 'COMMON.DERIV'
2059       include 'COMMON.GEO'
2060       include 'COMMON.LOCAL'
2061       include 'COMMON.INTERACT'
2062       include 'COMMON.IOUNITS'
2063       include 'COMMON.NAMES'
2064       include 'COMMON.TIME1'
2065       double precision xv,sigv,lowb,highb
2066 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2067 c First generate velocities in the eigenspace of the G matrix
2068 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2069 c      call flush(iout)
2070 c      write (iout,*) "RANDOM_VEL dimen",dimen
2071       xv=0.0d0
2072       ii=0
2073       do i=1,dimen
2074         do k=1,3
2075           ii=ii+1
2076           sigv=dsqrt((Rb*t_bath)/geigen(i))
2077           lowb=-5*sigv
2078           highb=5*sigv
2079           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2080 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2081 c     &      " d_t_work_new",d_t_work_new(ii)
2082         enddo
2083       enddo
2084       call flush(iout)
2085 c diagnostics
2086 c      Ek1=0.0d0
2087 c      ii=0
2088 c      do i=1,dimen
2089 c        do k=1,3
2090 c          ii=ii+1
2091 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2092 c        enddo
2093 c      enddo
2094 c      write (iout,*) "Ek from eigenvectors",Ek1
2095 c end diagnostics
2096 c Transform velocities to UNRES coordinate space
2097       do k=0,2
2098         do i=1,dimen
2099           ind=(i-1)*3+k+1
2100           d_t_work(ind)=0.0d0
2101           do j=1,dimen
2102             d_t_work(ind)=d_t_work(ind)
2103      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2104           enddo
2105 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2106 c          call flush(iout)
2107         enddo
2108       enddo
2109 c Transfer to the d_t vector
2110       do j=1,3
2111         d_t(j,0)=d_t_work(j)
2112       enddo
2113       ind=3
2114       do i=nnt,nct-1
2115         do j=1,3
2116           ind=ind+1
2117           d_t(j,i)=d_t_work(ind)
2118         enddo
2119       enddo
2120 c      do i=0,nres-1
2121 c        write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2122 c      enddo
2123 c      call flush(iout)
2124       do i=nnt,nct
2125         if (itype(i).ne.10) then
2126           do j=1,3
2127             ind=ind+1
2128             d_t(j,i+nres)=d_t_work(ind)
2129           enddo
2130         endif
2131       enddo
2132 c      call kinetic(EK)
2133 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2134 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2135 c      call flush(iout)
2136       return
2137       end
2138 #ifndef LANG0
2139 c-----------------------------------------------------------
2140       subroutine sd_verlet_p_setup
2141 c Sets up the parameters of stochastic Verlet algorithm
2142       implicit real*8 (a-h,o-z)
2143       include 'DIMENSIONS'
2144 #ifdef MPI
2145       include 'mpif.h'
2146 #endif
2147       include 'COMMON.CONTROL'
2148       include 'COMMON.VAR'
2149       include 'COMMON.MD'
2150 #ifndef LANG0
2151       include 'COMMON.LANGEVIN'
2152 #else
2153       include 'COMMON.LANGEVIN.lang0'
2154 #endif
2155       include 'COMMON.CHAIN'
2156       include 'COMMON.DERIV'
2157       include 'COMMON.GEO'
2158       include 'COMMON.LOCAL'
2159       include 'COMMON.INTERACT'
2160       include 'COMMON.IOUNITS'
2161       include 'COMMON.NAMES'
2162       include 'COMMON.TIME1'
2163       double precision emgdt(MAXRES6),
2164      & pterm,vterm,rho,rhoc,vsig,
2165      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2166      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2167      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2168       logical lprn /.false./
2169       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2170       double precision ktm
2171 #ifdef MPI
2172       tt0 = MPI_Wtime()
2173 #else
2174       tt0 = tcpu()
2175 #endif
2176 c
2177 c AL 8/17/04 Code adapted from tinker
2178 c
2179 c Get the frictional and random terms for stochastic dynamics in the
2180 c eigenspace of mass-scaled UNRES friction matrix
2181 c
2182       do i = 1, dimen
2183             gdt = fricgam(i) * d_time
2184 c
2185 c Stochastic dynamics reduces to simple MD for zero friction
2186 c
2187             if (gdt .le. zero) then
2188                pfric_vec(i) = 1.0d0
2189                vfric_vec(i) = d_time
2190                afric_vec(i) = 0.5d0 * d_time * d_time
2191                prand_vec(i) = 0.0d0
2192                vrand_vec1(i) = 0.0d0
2193                vrand_vec2(i) = 0.0d0
2194 c
2195 c Analytical expressions when friction coefficient is large
2196 c
2197             else
2198                if (gdt .ge. gdt_radius) then
2199                   egdt = dexp(-gdt)
2200                   pfric_vec(i) = egdt
2201                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2202                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2203                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2204                   vterm = 1.0d0 - egdt**2
2205                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2206 c
2207 c Use series expansions when friction coefficient is small
2208 c
2209                else
2210                   gdt2 = gdt * gdt
2211                   gdt3 = gdt * gdt2
2212                   gdt4 = gdt2 * gdt2
2213                   gdt5 = gdt2 * gdt3
2214                   gdt6 = gdt3 * gdt3
2215                   gdt7 = gdt3 * gdt4
2216                   gdt8 = gdt4 * gdt4
2217                   gdt9 = gdt4 * gdt5
2218                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2219      &                          - gdt5/120.0d0 + gdt6/720.0d0
2220      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2221      &                          - gdt9/362880.0d0) / fricgam(i)**2
2222                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2223                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2224                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2225      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2226      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2227      &                       + 127.0d0*gdt9/90720.0d0
2228                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2229      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2230      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2231      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2232                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2233      &                       - 17.0d0*gdt2/1280.0d0
2234      &                       + 17.0d0*gdt3/6144.0d0
2235      &                       + 40967.0d0*gdt4/34406400.0d0
2236      &                       - 57203.0d0*gdt5/275251200.0d0
2237      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2238                end if
2239 c
2240 c Compute the scaling factors of random terms for the nonzero friction case
2241 c
2242                ktm = 0.5d0*d_time/fricgam(i)
2243                psig = dsqrt(ktm*pterm) / fricgam(i)
2244                vsig = dsqrt(ktm*vterm)
2245                rhoc = dsqrt(1.0d0 - rho*rho)
2246                prand_vec(i) = psig
2247                vrand_vec1(i) = vsig * rho
2248                vrand_vec2(i) = vsig * rhoc
2249             end if
2250       end do
2251       if (lprn) then
2252       write (iout,*)
2253      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2254      &  " vrand_vec2"
2255       do i=1,dimen
2256         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2257      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2258       enddo
2259       endif
2260 c
2261 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2262 c
2263 #ifndef   LANG0
2264       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2265       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2266       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2267       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2268       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2269       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2270 #endif
2271 #ifdef MPI
2272       t_sdsetup=t_sdsetup+MPI_Wtime()
2273 #else
2274       t_sdsetup=t_sdsetup+tcpu()-tt0
2275 #endif
2276       return
2277       end
2278 c-------------------------------------------------------------
2279       subroutine eigtransf1(n,ndim,ab,d,c)
2280       implicit none
2281       integer n,ndim
2282       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2283       integer i,j,k
2284       do i=1,n
2285         do j=1,n
2286           c(i,j)=0.0d0
2287           do k=1,n
2288             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2289           enddo
2290         enddo
2291       enddo
2292       return
2293       end
2294 c-------------------------------------------------------------
2295       subroutine eigtransf(n,ndim,a,b,d,c)
2296       implicit none
2297       integer n,ndim
2298       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2299       integer i,j,k
2300       do i=1,n
2301         do j=1,n
2302           c(i,j)=0.0d0
2303           do k=1,n
2304             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2305           enddo
2306         enddo
2307       enddo
2308       return
2309       end
2310 c-------------------------------------------------------------
2311       subroutine sd_verlet1
2312 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2313       implicit real*8 (a-h,o-z)
2314       include 'DIMENSIONS'
2315       include 'COMMON.CONTROL'
2316       include 'COMMON.VAR'
2317       include 'COMMON.MD'
2318 #ifndef LANG0
2319       include 'COMMON.LANGEVIN'
2320 #else
2321       include 'COMMON.LANGEVIN.lang0'
2322 #endif
2323       include 'COMMON.CHAIN'
2324       include 'COMMON.DERIV'
2325       include 'COMMON.GEO'
2326       include 'COMMON.LOCAL'
2327       include 'COMMON.INTERACT'
2328       include 'COMMON.IOUNITS'
2329       include 'COMMON.NAMES'
2330       double precision stochforcvec(MAXRES6)
2331       common /stochcalc/ stochforcvec
2332       logical lprn /.false./
2333
2334 c      write (iout,*) "dc_old"
2335 c      do i=0,nres
2336 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2337 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2338 c      enddo
2339       do j=1,3
2340         dc_work(j)=dc_old(j,0)
2341         d_t_work(j)=d_t_old(j,0)
2342         d_a_work(j)=d_a_old(j,0)
2343       enddo
2344       ind=3
2345       do i=nnt,nct-1
2346         do j=1,3
2347           dc_work(ind+j)=dc_old(j,i)
2348           d_t_work(ind+j)=d_t_old(j,i)
2349           d_a_work(ind+j)=d_a_old(j,i)
2350         enddo
2351         ind=ind+3
2352       enddo
2353       do i=nnt,nct
2354         if (itype(i).ne.10) then
2355           do j=1,3
2356             dc_work(ind+j)=dc_old(j,i+nres)
2357             d_t_work(ind+j)=d_t_old(j,i+nres)
2358             d_a_work(ind+j)=d_a_old(j,i+nres)
2359           enddo
2360           ind=ind+3
2361         endif
2362       enddo
2363 #ifndef LANG0
2364       if (lprn) then
2365       write (iout,*)
2366      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2367      &  " vrand_mat2"
2368       do i=1,dimen
2369         do j=1,dimen
2370           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2371      &      vfric_mat(i,j),afric_mat(i,j),
2372      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2373         enddo
2374       enddo
2375       endif
2376       do i=1,dimen
2377         ddt1=0.0d0
2378         ddt2=0.0d0
2379         do j=1,dimen
2380           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2381      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2382           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2383           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2384         enddo
2385         d_t_work_new(i)=ddt1+0.5d0*ddt2
2386         d_t_work(i)=ddt1+ddt2
2387       enddo
2388 #endif
2389       do j=1,3
2390         dc(j,0)=dc_work(j)
2391         d_t(j,0)=d_t_work(j)
2392       enddo
2393       ind=3
2394       do i=nnt,nct-1
2395         do j=1,3
2396           dc(j,i)=dc_work(ind+j)
2397           d_t(j,i)=d_t_work(ind+j)
2398         enddo
2399         ind=ind+3
2400       enddo
2401       do i=nnt,nct
2402         if (itype(i).ne.10) then
2403           inres=i+nres
2404           do j=1,3
2405             dc(j,inres)=dc_work(ind+j)
2406             d_t(j,inres)=d_t_work(ind+j)
2407           enddo
2408           ind=ind+3
2409         endif
2410       enddo
2411       return
2412       end
2413 c--------------------------------------------------------------------------
2414       subroutine sd_verlet2
2415 c  Calculating the adjusted velocities for accelerations
2416       implicit real*8 (a-h,o-z)
2417       include 'DIMENSIONS'
2418       include 'COMMON.CONTROL'
2419       include 'COMMON.VAR'
2420       include 'COMMON.MD'
2421 #ifndef LANG0
2422       include 'COMMON.LANGEVIN'
2423 #else
2424       include 'COMMON.LANGEVIN.lang0'
2425 #endif
2426       include 'COMMON.CHAIN'
2427       include 'COMMON.DERIV'
2428       include 'COMMON.GEO'
2429       include 'COMMON.LOCAL'
2430       include 'COMMON.INTERACT'
2431       include 'COMMON.IOUNITS'
2432       include 'COMMON.NAMES'
2433       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2434       common /stochcalc/ stochforcvec
2435 c
2436 c Compute the stochastic forces which contribute to velocity change
2437 c
2438       call stochastic_force(stochforcvecV)
2439
2440 #ifndef LANG0
2441       do i=1,dimen
2442         ddt1=0.0d0
2443         ddt2=0.0d0
2444         do j=1,dimen
2445           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2446           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2447      &     vrand_mat2(i,j)*stochforcvecV(j)
2448         enddo
2449         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2450       enddo
2451 #endif
2452       do j=1,3
2453         d_t(j,0)=d_t_work(j)
2454       enddo
2455       ind=3
2456       do i=nnt,nct-1
2457         do j=1,3
2458           d_t(j,i)=d_t_work(ind+j)
2459         enddo
2460         ind=ind+3
2461       enddo
2462       do i=nnt,nct
2463         if (itype(i).ne.10) then
2464           inres=i+nres
2465           do j=1,3
2466             d_t(j,inres)=d_t_work(ind+j)
2467           enddo
2468           ind=ind+3
2469         endif
2470       enddo
2471       return
2472       end
2473 c-----------------------------------------------------------
2474       subroutine sd_verlet_ciccotti_setup
2475 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2476 c version
2477       implicit real*8 (a-h,o-z)
2478       include 'DIMENSIONS'
2479 #ifdef MPI
2480       include 'mpif.h'
2481 #endif
2482       include 'COMMON.CONTROL'
2483       include 'COMMON.VAR'
2484       include 'COMMON.MD'
2485 #ifndef LANG0
2486       include 'COMMON.LANGEVIN'
2487 #else
2488       include 'COMMON.LANGEVIN.lang0'
2489 #endif
2490       include 'COMMON.CHAIN'
2491       include 'COMMON.DERIV'
2492       include 'COMMON.GEO'
2493       include 'COMMON.LOCAL'
2494       include 'COMMON.INTERACT'
2495       include 'COMMON.IOUNITS'
2496       include 'COMMON.NAMES'
2497       include 'COMMON.TIME1'
2498       double precision emgdt(MAXRES6),
2499      & pterm,vterm,rho,rhoc,vsig,
2500      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2501      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2502      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2503       logical lprn /.false./
2504       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2505       double precision ktm
2506 #ifdef MPI
2507       tt0 = MPI_Wtime()
2508 #else
2509       tt0 = tcpu()
2510 #endif
2511 c
2512 c AL 8/17/04 Code adapted from tinker
2513 c
2514 c Get the frictional and random terms for stochastic dynamics in the
2515 c eigenspace of mass-scaled UNRES friction matrix
2516 c
2517       do i = 1, dimen
2518             write (iout,*) "i",i," fricgam",fricgam(i)
2519             gdt = fricgam(i) * d_time
2520 c
2521 c Stochastic dynamics reduces to simple MD for zero friction
2522 c
2523             if (gdt .le. zero) then
2524                pfric_vec(i) = 1.0d0
2525                vfric_vec(i) = d_time
2526                afric_vec(i) = 0.5d0*d_time*d_time
2527                prand_vec(i) = afric_vec(i)
2528                vrand_vec2(i) = vfric_vec(i)
2529 c
2530 c Analytical expressions when friction coefficient is large
2531 c
2532             else
2533                egdt = dexp(-gdt)
2534                pfric_vec(i) = egdt
2535                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2536                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2537                prand_vec(i) = afric_vec(i)
2538                vrand_vec2(i) = vfric_vec(i)
2539 c
2540 c Compute the scaling factors of random terms for the nonzero friction case
2541 c
2542 c               ktm = 0.5d0*d_time/fricgam(i)
2543 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2544 c               vsig = dsqrt(ktm*vterm)
2545 c               prand_vec(i) = psig*afric_vec(i)
2546 c               vrand_vec2(i) = vsig*vfric_vec(i)
2547             end if
2548       end do
2549       if (lprn) then
2550       write (iout,*)
2551      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2552      &  " vrand_vec2"
2553       do i=1,dimen
2554         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2555      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2556       enddo
2557       endif
2558 c
2559 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2560 c
2561       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2562       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2563       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2564       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2565       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2566 #ifdef MPI
2567       t_sdsetup=t_sdsetup+MPI_Wtime()
2568 #else
2569       t_sdsetup=t_sdsetup+tcpu()-tt0
2570 #endif
2571       return
2572       end
2573 c-------------------------------------------------------------
2574       subroutine sd_verlet1_ciccotti
2575 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2576       implicit real*8 (a-h,o-z)
2577       include 'DIMENSIONS'
2578 #ifdef MPI
2579       include 'mpif.h'
2580 #endif
2581       include 'COMMON.CONTROL'
2582       include 'COMMON.VAR'
2583       include 'COMMON.MD'
2584 #ifndef LANG0
2585       include 'COMMON.LANGEVIN'
2586 #else
2587       include 'COMMON.LANGEVIN.lang0'
2588 #endif
2589       include 'COMMON.CHAIN'
2590       include 'COMMON.DERIV'
2591       include 'COMMON.GEO'
2592       include 'COMMON.LOCAL'
2593       include 'COMMON.INTERACT'
2594       include 'COMMON.IOUNITS'
2595       include 'COMMON.NAMES'
2596       double precision stochforcvec(MAXRES6)
2597       common /stochcalc/ stochforcvec
2598       logical lprn /.false./
2599
2600 c      write (iout,*) "dc_old"
2601 c      do i=0,nres
2602 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2603 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2604 c      enddo
2605       do j=1,3
2606         dc_work(j)=dc_old(j,0)
2607         d_t_work(j)=d_t_old(j,0)
2608         d_a_work(j)=d_a_old(j,0)
2609       enddo
2610       ind=3
2611       do i=nnt,nct-1
2612         do j=1,3
2613           dc_work(ind+j)=dc_old(j,i)
2614           d_t_work(ind+j)=d_t_old(j,i)
2615           d_a_work(ind+j)=d_a_old(j,i)
2616         enddo
2617         ind=ind+3
2618       enddo
2619       do i=nnt,nct
2620         if (itype(i).ne.10) then
2621           do j=1,3
2622             dc_work(ind+j)=dc_old(j,i+nres)
2623             d_t_work(ind+j)=d_t_old(j,i+nres)
2624             d_a_work(ind+j)=d_a_old(j,i+nres)
2625           enddo
2626           ind=ind+3
2627         endif
2628       enddo
2629
2630 #ifndef LANG0
2631       if (lprn) then
2632       write (iout,*)
2633      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2634      &  " vrand_mat2"
2635       do i=1,dimen
2636         do j=1,dimen
2637           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2638      &      vfric_mat(i,j),afric_mat(i,j),
2639      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2640         enddo
2641       enddo
2642       endif
2643       do i=1,dimen
2644         ddt1=0.0d0
2645         ddt2=0.0d0
2646         do j=1,dimen
2647           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2648      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2649           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2650           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2651         enddo
2652         d_t_work_new(i)=ddt1+0.5d0*ddt2
2653         d_t_work(i)=ddt1+ddt2
2654       enddo
2655 #endif
2656       do j=1,3
2657         dc(j,0)=dc_work(j)
2658         d_t(j,0)=d_t_work(j)
2659       enddo
2660       ind=3
2661       do i=nnt,nct-1
2662         do j=1,3
2663           dc(j,i)=dc_work(ind+j)
2664           d_t(j,i)=d_t_work(ind+j)
2665         enddo
2666         ind=ind+3
2667       enddo
2668       do i=nnt,nct
2669         if (itype(i).ne.10) then
2670           inres=i+nres
2671           do j=1,3
2672             dc(j,inres)=dc_work(ind+j)
2673             d_t(j,inres)=d_t_work(ind+j)
2674           enddo
2675           ind=ind+3
2676         endif
2677       enddo
2678       return
2679       end
2680 c--------------------------------------------------------------------------
2681       subroutine sd_verlet2_ciccotti
2682 c  Calculating the adjusted velocities for accelerations
2683       implicit real*8 (a-h,o-z)
2684       include 'DIMENSIONS'
2685       include 'COMMON.CONTROL'
2686       include 'COMMON.VAR'
2687       include 'COMMON.MD'
2688 #ifndef LANG0
2689       include 'COMMON.LANGEVIN'
2690 #else
2691       include 'COMMON.LANGEVIN.lang0'
2692 #endif
2693       include 'COMMON.CHAIN'
2694       include 'COMMON.DERIV'
2695       include 'COMMON.GEO'
2696       include 'COMMON.LOCAL'
2697       include 'COMMON.INTERACT'
2698       include 'COMMON.IOUNITS'
2699       include 'COMMON.NAMES'
2700       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2701       common /stochcalc/ stochforcvec
2702 c
2703 c Compute the stochastic forces which contribute to velocity change
2704 c
2705       call stochastic_force(stochforcvecV)
2706 #ifndef LANG0
2707       do i=1,dimen
2708         ddt1=0.0d0
2709         ddt2=0.0d0
2710         do j=1,dimen
2711
2712           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2713 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2714           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2715         enddo
2716         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2717       enddo
2718 #endif
2719       do j=1,3
2720         d_t(j,0)=d_t_work(j)
2721       enddo
2722       ind=3
2723       do i=nnt,nct-1
2724         do j=1,3
2725           d_t(j,i)=d_t_work(ind+j)
2726         enddo
2727         ind=ind+3
2728       enddo
2729       do i=nnt,nct
2730         if (itype(i).ne.10) then
2731           inres=i+nres
2732           do j=1,3
2733             d_t(j,inres)=d_t_work(ind+j)
2734           enddo
2735           ind=ind+3
2736         endif
2737       enddo
2738       return
2739       end
2740 #endif
2741 c------------------------------------------------------
2742       double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2743       implicit none
2744       double precision ek,s,e,pi,Q,t_bath,Rb
2745       integer dimenl
2746       Rb=0.001986d0
2747       HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2748 c      print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2749 c     &      pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2750       return
2751       end
2752 c-----------------------------------------------------------------
2753       double precision function HNose_nh(eki,e)
2754       implicit real*8 (a-h,o-z)
2755       include 'DIMENSIONS'
2756       include 'COMMON.MD'
2757       HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2758       do i=2,nnos
2759         HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2760       enddo
2761 c      write(4,'(5e15.5)')
2762 c     &       vlogs(1),xlogs(1),HNose,eki,e
2763       return
2764       end
2765 c-----------------------------------------------------------------
2766       SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2767       implicit real*8 (a-h,o-z)
2768       include 'DIMENSIONS'
2769       include 'COMMON.MD'
2770       double precision akin,gnkt,dt,aa,gkt,scale
2771       double precision wdti(maxyosh),wdti2(maxyosh),
2772      &                 wdti4(maxyosh),wdti8(maxyosh)
2773       integer i,iresn,iyosh,inos,nnos1
2774
2775       dt=d_time
2776       nnos1=nnos+1
2777       GKT = Rb*t_bath
2778       GNKT = dimen3*GKT
2779       akin=akin*2
2780
2781
2782 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2783 C INTEGRATION FROM t=0 TO t=DT/2
2784 C GET THE TOTAL KINETIC ENERGY
2785       SCALE = 1.D0
2786 c      CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2787 C UPDATE THE FORCES
2788       GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2789 C START THE MULTIPLE TIME STEP PROCEDURE
2790       DO IRESN = 1,NRESN
2791        DO IYOSH = 1,NYOSH
2792 C UPDATE THE THERMOSTAT VELOCITIES
2793         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2794         DO INOS = 1,NNOS-1
2795          AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2796          VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2797      &          + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2798         ENDDO
2799 C UPDATE THE PARTICLE VELOCITIES
2800         AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2801         SCALE = SCALE*AA
2802 C UPDATE THE FORCES
2803         GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2804 C UPDATE THE THERMOSTAT POSITIONS
2805         DO INOS = 1,NNOS
2806          XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2807         ENDDO
2808 C UPDATE THE THERMOSTAT VELOCITIES
2809         DO INOS = 1,NNOS-1
2810          AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2811          VLOGS(INOS) = VLOGS(INOS)*AA*AA
2812      &      + WDTI4(IYOSH)*GLOGS(INOS)*AA
2813          GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2814      &      -GKT)/QMASS(INOS+1)
2815         ENDDO
2816         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2817        ENDDO
2818       ENDDO
2819 C UPDATE THE PARTICLE VELOCITIES
2820 c outside of this subroutine
2821 c      DO I = 1,N
2822 c       VX(I) = VX(I)*SCALE
2823 c       VY(I) = VY(I)*SCALE
2824 c       VZ(I) = VZ(I)*SCALE
2825 c      ENDDO
2826       RETURN
2827       END
2828 c-----------------------------------------------------------------
2829       subroutine tnp1_respa_i_step1
2830 c Applying Nose-Poincare algorithm - step 1 to coordinates
2831 c JPSJ 70 75 (2001) S. Nose
2832 c
2833 c d_t is not updated here
2834 c
2835       implicit real*8 (a-h,o-z)
2836       include 'DIMENSIONS'
2837       include 'COMMON.CONTROL'
2838       include 'COMMON.VAR'
2839       include 'COMMON.MD'
2840       include 'COMMON.CHAIN'
2841       include 'COMMON.DERIV'
2842       include 'COMMON.GEO'
2843       include 'COMMON.LOCAL'
2844       include 'COMMON.INTERACT'
2845       include 'COMMON.IOUNITS'
2846       include 'COMMON.NAMES'
2847       double precision adt,adt2,tmp
2848
2849       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2850       s12_np=s_np*tmp**2
2851       pistar=pi_np/tmp
2852       s12_dt=d_time/s12_np
2853       d_time_s12=d_time*0.5*s12_np
2854
2855       do j=1,3
2856         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2857         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2858       enddo
2859       do i=nnt,nct-1
2860         do j=1,3
2861           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2862           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2863         enddo
2864       enddo
2865       do i=nnt,nct
2866         if (itype(i).ne.10) then
2867           inres=i+nres
2868           do j=1,3
2869            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2870            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2871           enddo
2872         endif
2873       enddo
2874       return
2875       end
2876 c---------------------------------------------------------------------
2877       subroutine tnp1_respa_i_step2
2878 c  Step 2 of the velocity Verlet algorithm: update velocities
2879       implicit real*8 (a-h,o-z)
2880       include 'DIMENSIONS'
2881       include 'COMMON.CONTROL'
2882       include 'COMMON.VAR'
2883       include 'COMMON.MD'
2884       include 'COMMON.CHAIN'
2885       include 'COMMON.DERIV'
2886       include 'COMMON.GEO'
2887       include 'COMMON.LOCAL'
2888       include 'COMMON.INTERACT'
2889       include 'COMMON.IOUNITS'
2890       include 'COMMON.NAMES'
2891
2892       double precision d_time_s12
2893
2894       do i=0,2*nres
2895        do j=1,3
2896         d_t(j,i)=d_t_new(j,i)
2897        enddo
2898       enddo
2899
2900       call kinetic(EK)
2901       EK=EK/s12_np**2
2902
2903       d_time_s12=0.5d0*s12_np*d_time
2904
2905       do j=1,3
2906         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2907       enddo
2908       do i=nnt,nct-1
2909         do j=1,3
2910           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2911         enddo
2912       enddo
2913       do i=nnt,nct
2914         if (itype(i).ne.10) then
2915           inres=i+nres
2916           do j=1,3
2917             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2918           enddo
2919         endif
2920       enddo
2921
2922       pistar=pistar+(EK-0.5*(E_old+potE)
2923      &       -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2924       tmp=1+pistar/(2*Q_np)*0.5*d_time
2925       s_np=s12_np*tmp**2
2926       pi_np=pistar/tmp
2927
2928       return
2929       end
2930 c-------------------------------------------------------
2931
2932       subroutine tnp1_step1
2933 c Applying Nose-Poincare algorithm - step 1 to coordinates
2934 c JPSJ 70 75 (2001) S. Nose
2935 c
2936 c d_t is not updated here
2937 c
2938       implicit real*8 (a-h,o-z)
2939       include 'DIMENSIONS'
2940       include 'COMMON.CONTROL'
2941       include 'COMMON.VAR'
2942       include 'COMMON.MD'
2943       include 'COMMON.CHAIN'
2944       include 'COMMON.DERIV'
2945       include 'COMMON.GEO'
2946       include 'COMMON.LOCAL'
2947       include 'COMMON.INTERACT'
2948       include 'COMMON.IOUNITS'
2949       include 'COMMON.NAMES'
2950       double precision adt,adt2,tmp
2951
2952       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2953       s12_np=s_np*tmp**2
2954       pistar=pi_np/tmp
2955       s12_dt=d_time/s12_np
2956       d_time_s12=d_time*0.5*s12_np
2957
2958       do j=1,3
2959         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2960         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2961       enddo
2962       do i=nnt,nct-1
2963         do j=1,3
2964           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2965           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2966         enddo
2967       enddo
2968       do i=nnt,nct
2969         if (itype(i).ne.10) then
2970           inres=i+nres
2971           do j=1,3
2972            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2973            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2974           enddo
2975         endif
2976       enddo
2977       return
2978       end
2979 c---------------------------------------------------------------------
2980       subroutine tnp1_step2
2981 c  Step 2 of the velocity Verlet algorithm: update velocities
2982       implicit real*8 (a-h,o-z)
2983       include 'DIMENSIONS'
2984       include 'COMMON.CONTROL'
2985       include 'COMMON.VAR'
2986       include 'COMMON.MD'
2987       include 'COMMON.CHAIN'
2988       include 'COMMON.DERIV'
2989       include 'COMMON.GEO'
2990       include 'COMMON.LOCAL'
2991       include 'COMMON.INTERACT'
2992       include 'COMMON.IOUNITS'
2993       include 'COMMON.NAMES'
2994
2995       double precision d_time_s12
2996
2997       do i=0,2*nres
2998        do j=1,3
2999         d_t(j,i)=d_t_new(j,i)
3000        enddo
3001       enddo
3002
3003       call kinetic(EK)
3004       EK=EK/s12_np**2
3005
3006       d_time_s12=0.5d0*s12_np*d_time
3007
3008       do j=1,3
3009         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3010       enddo
3011       do i=nnt,nct-1
3012         do j=1,3
3013           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3014         enddo
3015       enddo
3016       do i=nnt,nct
3017         if (itype(i).ne.10) then
3018           inres=i+nres
3019           do j=1,3
3020             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3021           enddo
3022         endif
3023       enddo
3024
3025 cd      write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3026       pistar=pistar+(EK-0.5*(E_old+potE)
3027      &       -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3028       tmp=1+pistar/(2*Q_np)*0.5*d_time
3029       s_np=s12_np*tmp**2
3030       pi_np=pistar/tmp
3031
3032       return
3033       end
3034
3035 c-----------------------------------------------------------------
3036       subroutine tnp_respa_i_step1
3037 c Applying Nose-Poincare algorithm - step 1 to coordinates
3038 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3039 c
3040 c d_t is not updated here, it is destroyed
3041 c
3042       implicit real*8 (a-h,o-z)
3043       include 'DIMENSIONS'
3044       include 'COMMON.CONTROL'
3045       include 'COMMON.VAR'
3046       include 'COMMON.MD'
3047       include 'COMMON.CHAIN'
3048       include 'COMMON.DERIV'
3049       include 'COMMON.GEO'
3050       include 'COMMON.LOCAL'
3051       include 'COMMON.INTERACT'
3052       include 'COMMON.IOUNITS'
3053       include 'COMMON.NAMES'
3054       double precision C_np,d_time_s,tmp,d_time_ss
3055
3056       d_time_s=d_time*0.5*s_np
3057 ct2      d_time_s=d_time*0.5*s12_np
3058
3059       do j=1,3
3060         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3061       enddo
3062       do i=nnt,nct-1
3063         do j=1,3
3064           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3065         enddo
3066       enddo
3067       do i=nnt,nct
3068         if (itype(i).ne.10) then
3069           inres=i+nres
3070           do j=1,3
3071            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3072           enddo
3073         endif
3074       enddo
3075
3076       do i=0,2*nres
3077        do j=1,3
3078         d_t(j,i)=d_t_new(j,i)
3079        enddo
3080       enddo
3081
3082       call kinetic(EK)
3083       EK=EK/s_np**2
3084
3085       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3086      &                     -pi_np
3087
3088       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3089       tmp=0.5*d_time*pistar/Q_np
3090       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3091
3092       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3093 ct2      d_time_ss=d_time/s12_np
3094 c      d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3095
3096       do j=1,3
3097         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3098       enddo
3099       do i=nnt,nct-1
3100         do j=1,3
3101           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3102         enddo
3103       enddo
3104       do i=nnt,nct
3105         if (itype(i).ne.10) then
3106           inres=i+nres
3107           do j=1,3
3108            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3109           enddo
3110         endif
3111       enddo
3112
3113       return
3114       end
3115 c---------------------------------------------------------------------
3116
3117       subroutine tnp_respa_i_step2
3118 c  Step 2 of the velocity Verlet algorithm: update velocities
3119       implicit real*8 (a-h,o-z)
3120       include 'DIMENSIONS'
3121       include 'COMMON.CONTROL'
3122       include 'COMMON.VAR'
3123       include 'COMMON.MD'
3124       include 'COMMON.CHAIN'
3125       include 'COMMON.DERIV'
3126       include 'COMMON.GEO'
3127       include 'COMMON.LOCAL'
3128       include 'COMMON.INTERACT'
3129       include 'COMMON.IOUNITS'
3130       include 'COMMON.NAMES'
3131
3132       double precision d_time_s
3133
3134       EK=EK*(s_np/s12_np)**2
3135       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3136       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3137      &                              -HNose1+Csplit)
3138
3139 cr      print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3140       d_time_s=d_time*0.5*s12_np
3141 c      d_time_s=d_time*0.5*s_np
3142
3143       do j=1,3
3144         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3145       enddo
3146       do i=nnt,nct-1
3147         do j=1,3
3148           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3149         enddo
3150       enddo
3151       do i=nnt,nct
3152         if (itype(i).ne.10) then
3153           inres=i+nres
3154           do j=1,3
3155             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3156           enddo
3157         endif
3158       enddo
3159
3160       s_np=s12_np
3161
3162       return
3163       end
3164 c-----------------------------------------------------------------
3165       subroutine tnp_respa_step1
3166 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3167 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3168 c
3169 c d_t is not updated here, it is destroyed
3170 c
3171       implicit real*8 (a-h,o-z)
3172       include 'DIMENSIONS'
3173       include 'COMMON.CONTROL'
3174       include 'COMMON.VAR'
3175       include 'COMMON.MD'
3176       include 'COMMON.CHAIN'
3177       include 'COMMON.DERIV'
3178       include 'COMMON.GEO'
3179       include 'COMMON.LOCAL'
3180       include 'COMMON.INTERACT'
3181       include 'COMMON.IOUNITS'
3182       include 'COMMON.NAMES'
3183       double precision C_np,d_time_s,tmp,d_time_ss
3184       double precision energia(0:n_ene)
3185
3186       d_time_s=d_time*0.5*s_np
3187
3188       do j=1,3
3189         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3190       enddo
3191       do i=nnt,nct-1
3192         do j=1,3
3193           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3194         enddo
3195       enddo
3196       do i=nnt,nct
3197         if (itype(i).ne.10) then
3198           inres=i+nres
3199           do j=1,3
3200            d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3201           enddo
3202         endif
3203       enddo
3204
3205
3206 c      C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3207 c     &                     -pi_np
3208 c
3209 c      pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3210 c      tmp=0.5*d_time*pistar/Q_np
3211 c      s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3212 c      write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3213
3214 ct1      pi_np=pistar
3215 c      sold_np=s_np
3216 c      s_np=s12_np
3217
3218 c-------------------------------------
3219 c test of reviewer's comment
3220        pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3221 cr       print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3222 c-------------------------------------
3223
3224       return
3225       end
3226 c---------------------------------------------------------------------
3227       subroutine tnp_respa_step2
3228 c  Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3229       implicit real*8 (a-h,o-z)
3230       include 'DIMENSIONS'
3231       include 'COMMON.CONTROL'
3232       include 'COMMON.VAR'
3233       include 'COMMON.MD'
3234       include 'COMMON.CHAIN'
3235       include 'COMMON.DERIV'
3236       include 'COMMON.GEO'
3237       include 'COMMON.LOCAL'
3238       include 'COMMON.INTERACT'
3239       include 'COMMON.IOUNITS'
3240       include 'COMMON.NAMES'
3241
3242       double precision d_time_s
3243
3244 ct1      s12_np=s_np
3245 ct2      pistar=pi_np
3246
3247 ct      call kinetic(EK)
3248 ct      HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3249 ct      pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3250 ct     &                              -0.5*d_time*(HNose1-H0)
3251
3252 c-------------------------------------
3253 c test of reviewer's comment
3254       pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3255 cr      print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3256 c-------------------------------------
3257       d_time_s=d_time*0.5*s_np
3258
3259       do j=1,3
3260         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3261       enddo
3262       do i=nnt,nct-1
3263         do j=1,3
3264           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3265         enddo
3266       enddo
3267       do i=nnt,nct
3268         if (itype(i).ne.10) then
3269           inres=i+nres
3270           do j=1,3
3271             d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3272           enddo
3273         endif
3274       enddo
3275
3276 cd      s_np=s12_np
3277
3278       return
3279       end
3280 c---------------------------------------------------------------------
3281       subroutine tnp_step1
3282 c Applying Nose-Poincare algorithm - step 1 to coordinates
3283 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3284 c
3285 c d_t is not updated here, it is destroyed
3286 c
3287       implicit real*8 (a-h,o-z)
3288       include 'DIMENSIONS'
3289       include 'COMMON.CONTROL'
3290       include 'COMMON.VAR'
3291       include 'COMMON.MD'
3292       include 'COMMON.CHAIN'
3293       include 'COMMON.DERIV'
3294       include 'COMMON.GEO'
3295       include 'COMMON.LOCAL'
3296       include 'COMMON.INTERACT'
3297       include 'COMMON.IOUNITS'
3298       include 'COMMON.NAMES'
3299       double precision C_np,d_time_s,tmp,d_time_ss
3300
3301       d_time_s=d_time*0.5*s_np
3302
3303       do j=1,3
3304         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3305       enddo
3306       do i=nnt,nct-1
3307         do j=1,3
3308           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3309         enddo
3310       enddo
3311       do i=nnt,nct
3312         if (itype(i).ne.10) then
3313           inres=i+nres
3314           do j=1,3
3315            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3316           enddo
3317         endif
3318       enddo
3319
3320       do i=0,2*nres
3321        do j=1,3
3322         d_t(j,i)=d_t_new(j,i)
3323        enddo
3324       enddo
3325
3326       call kinetic(EK)
3327       EK=EK/s_np**2
3328
3329       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3330      &                     -pi_np
3331
3332       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3333       tmp=0.5*d_time*pistar/Q_np
3334       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3335 c      write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3336
3337       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3338
3339       do j=1,3
3340         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3341       enddo
3342       do i=nnt,nct-1
3343         do j=1,3
3344           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3345         enddo
3346       enddo
3347       do i=nnt,nct
3348         if (itype(i).ne.10) then
3349           inres=i+nres
3350           do j=1,3
3351            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3352           enddo
3353         endif
3354       enddo
3355
3356       return
3357       end
3358 c-----------------------------------------------------------------
3359       subroutine tnp_step2
3360 c  Step 2 of the velocity Verlet algorithm: update velocities
3361       implicit real*8 (a-h,o-z)
3362       include 'DIMENSIONS'
3363       include 'COMMON.CONTROL'
3364       include 'COMMON.VAR'
3365       include 'COMMON.MD'
3366       include 'COMMON.CHAIN'
3367       include 'COMMON.DERIV'
3368       include 'COMMON.GEO'
3369       include 'COMMON.LOCAL'
3370       include 'COMMON.INTERACT'
3371       include 'COMMON.IOUNITS'
3372       include 'COMMON.NAMES'
3373
3374       double precision d_time_s
3375
3376       EK=EK*(s_np/s12_np)**2
3377       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3378       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3379      &                              -0.5*d_time*(HNose1-H0)
3380
3381 cd      write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3382       d_time_s=d_time*0.5*s12_np
3383
3384       do j=1,3
3385         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3386       enddo
3387       do i=nnt,nct-1
3388         do j=1,3
3389           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3390         enddo
3391       enddo
3392       do i=nnt,nct
3393         if (itype(i).ne.10) then
3394           inres=i+nres
3395           do j=1,3
3396             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3397           enddo
3398         endif
3399       enddo
3400
3401       s_np=s12_np
3402
3403       return
3404       end
3405
3406       subroutine hmc_test(itime)
3407       implicit real*8 (a-h,o-z)
3408       include 'DIMENSIONS'
3409       include 'COMMON.CONTROL'
3410       include 'COMMON.MD'
3411       include 'COMMON.CHAIN'
3412
3413            hmc_acc=hmc_acc+1
3414            delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
3415            if (delta .lt. -50.0d0) then
3416                 delta=0.0d0
3417            else
3418                 delta=dexp(delta)
3419            endif
3420            xxx=ran_number(0.0d0,1.0d0)
3421
3422            if (me.eq.king .or. .not. out1file)
3423      &       write(iout,'(a8,i5,6f10.4)')
3424      &        'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
3425
3426            if (delta .le. xxx) then
3427             do i=1,2*nres
3428              do j=1,3
3429               dc(j,i)=dc_hmc(j,i)
3430              enddo
3431             enddo
3432             itime=itime-hmc
3433             totT=totThmc
3434            else
3435             if (me.eq.king .or. .not. out1file)
3436      &       write(iout,*) 'HMC accepting new'
3437             totThmc=totT
3438             do i=1,2*nres
3439              do j=1,3
3440               dc_hmc(j,i)=dc(j,i)
3441              enddo
3442             enddo
3443            endif
3444
3445            call chainbuild_cart
3446            call random_vel
3447            do i=0,2*nres
3448             do j=1,3
3449               d_t_old(j,i)=d_t(j,i)
3450             enddo
3451            enddo
3452            call kinetic(EK)
3453            kinetic_T=2.0d0/(dimen3*Rb)*EK
3454            call etotal(potEcomp)
3455            potE=potEcomp(0)
3456            hmc_etot=potE+EK
3457            if (me.eq.king .or. .not. out1file)
3458      &      write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK
3459
3460
3461       return
3462       end