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(itype(i))*dsqrt(gamsc(itype(i)))
1685         enddo
1686       endif
1687 c Open the pdb file for snapshotshots
1688 #ifdef MPI
1689       if(mdpdb) then
1690         if (ilen(tmpdir).gt.0)
1691      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1692      &      liczba(:ilen(liczba))//".pdb")
1693         open(ipdb,
1694      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1695      &  //".pdb")
1696       else
1697 #ifdef NOXDR
1698         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1699      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1700      &      liczba(:ilen(liczba))//".x")
1701         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1702      &  //".x"
1703 #else
1704         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1705      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1706      &      liczba(:ilen(liczba))//".cx")
1707         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1708      &  //".cx"
1709 #endif
1710       endif
1711 #else
1712       if(mdpdb) then
1713          if (ilen(tmpdir).gt.0)
1714      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1715          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1716       else
1717          if (ilen(tmpdir).gt.0)
1718      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1719          cartname=prefix(:ilen(prefix))//"_MD.cx"
1720       endif
1721 #endif
1722       if (usampl) then
1723         write (qstr,'(256(1h ))')
1724         ipos=1
1725         do i=1,nfrag
1726           iq = qinfrag(i,iset)*10
1727           iw = wfrag(i,iset)/100
1728           if (iw.gt.0) then
1729             if(me.eq.king.or..not.out1file)
1730      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1731             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1732             ipos=ipos+7
1733           endif
1734         enddo
1735         do i=1,npair
1736           iq = qinpair(i,iset)*10
1737           iw = wpair(i,iset)/100
1738           if (iw.gt.0) then
1739             if(me.eq.king.or..not.out1file)
1740      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1741             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1742             ipos=ipos+7
1743           endif
1744         enddo
1745 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1746 #ifdef NOXDR
1747 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1748 #else
1749 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1750 #endif
1751 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1752       endif
1753       icg=1
1754       if (rest) then
1755        if (restart1file) then
1756          if (me.eq.king)
1757      &     inquire(file=mremd_rst_name,exist=file_exist)
1758            write (*,*) me," Before broadcast: file_exist",file_exist
1759 #ifdef MPI
1760          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1761      &          IERR)
1762          write (*,*) me," After broadcast: file_exist",file_exist
1763 #endif
1764 c        inquire(file=mremd_rst_name,exist=file_exist)
1765         if(me.eq.king.or..not.out1file)
1766      &   write(iout,*) "Initial state read by master and distributed"
1767        else
1768          if (ilen(tmpdir).gt.0)
1769      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1770      &      //liczba(:ilen(liczba))//'.rst')
1771         inquire(file=rest2name,exist=file_exist)
1772        endif
1773        if(file_exist) then
1774          if(.not.restart1file) then
1775            if(me.eq.king.or..not.out1file)
1776      &      write(iout,*) "Initial state will be read from file ",
1777      &      rest2name(:ilen(rest2name))
1778            call readrst
1779          endif
1780          call rescale_weights(t_bath)
1781        else
1782         if(me.eq.king.or..not.out1file)then
1783          if (restart1file) then
1784           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1785      &       " does not exist"
1786          else
1787           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1788      &       " does not exist"
1789          endif
1790          write(iout,*) "Initial velocities randomly generated"
1791         endif
1792         call random_vel
1793         totT=0.0d0
1794        endif
1795       else
1796 c Generate initial velocities
1797         if(me.eq.king.or..not.out1file)
1798      &   write(iout,*) "Initial velocities randomly generated"
1799         call random_vel
1800         totT=0.0d0
1801       endif
1802 c      rest2name = prefix(:ilen(prefix))//'.rst'
1803       if(me.eq.king.or..not.out1file)then
1804        write (iout,*) "Initial velocities"
1805        do i=0,nres
1806          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1807      &   (d_t(j,i+nres),j=1,3)
1808        enddo
1809 c  Zeroing the total angular momentum of the system
1810        write(iout,*) "Calling the zero-angular
1811      & momentum subroutine"
1812       endif
1813       call inertia_tensor
1814 c  Getting the potential energy and forces and velocities and accelerations
1815       call vcm_vel(vcm)
1816 c      write (iout,*) "velocity of the center of the mass:"
1817 c      write (iout,*) (vcm(j),j=1,3)
1818       do j=1,3
1819         d_t(j,0)=d_t(j,0)-vcm(j)
1820       enddo
1821 c Removing the velocity of the center of mass
1822       call vcm_vel(vcm)
1823       if(me.eq.king.or..not.out1file)then
1824        write (iout,*) "vcm right after adjustment:"
1825        write (iout,*) (vcm(j),j=1,3)
1826       endif
1827       if (.not.rest) then
1828          call chainbuild
1829          if(iranconf.ne.0) then
1830           if (overlapsc) then
1831            print *, 'Calling OVERLAP_SC'
1832            call overlap_sc(fail)
1833           endif
1834
1835           if (searchsc) then
1836            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1837            print *,'SC_move',nft_sc,etot
1838            if(me.eq.king.or..not.out1file)
1839      &      write(iout,*) 'SC_move',nft_sc,etot
1840           endif
1841
1842           if(dccart)then
1843            print *, 'Calling MINIM_DC'
1844            call minim_dc(etot,iretcode,nfun)
1845           else
1846            call geom_to_var(nvar,varia)
1847            print *,'Calling MINIMIZE.'
1848            call minimize(etot,varia,iretcode,nfun)
1849            call var_to_geom(nvar,varia)
1850           endif
1851           if(me.eq.king.or..not.out1file)
1852      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1853          endif
1854       endif
1855       call chainbuild_cart
1856       call kinetic(EK)
1857       if (tbf) then
1858         call verlet_bath(EK)
1859       endif
1860       kinetic_T=2.0d0/(dimen3*Rb)*EK
1861       if(me.eq.king.or..not.out1file)then
1862        call cartprint
1863        call intout
1864       endif
1865 #ifdef MPI
1866       tt0=MPI_Wtime()
1867 #else
1868       tt0=tcpu()
1869 #endif
1870       call zerograd
1871       call etotal(potEcomp)
1872 #ifdef TIMING_ENE
1873 #ifdef MPI
1874       t_etotal=t_etotal+MPI_Wtime()-tt0
1875 #else
1876       t_etotal=t_etotal+tcpu()-tt0
1877 #endif
1878 #endif
1879       potE=potEcomp(0)
1880
1881       if(tnp .or. tnp1) then
1882        s_np=1.0
1883        pi_np=0.0
1884        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
1885        H0=Hnose1
1886        write(iout,*) 'H0= ',H0
1887       endif
1888
1889       if(tnh) then
1890        HNose1=Hnose_nh(EK,potE)
1891        H0=HNose1
1892        write (iout,*) 'H0= ',H0
1893       endif
1894
1895       if (hmc.gt.0) then
1896          hmc_acc=0
1897          hmc_etot=potE+EK
1898           if(me.eq.king.or..not.out1file)
1899      &       write(iout,*) 'HMC',hmc_etot,potE,EK
1900          do i=1,2*nres
1901            do j=1,3
1902             dc_hmc(j,i)=dc(j,i)
1903            enddo
1904          enddo
1905       endif
1906
1907       call cartgrad
1908       call lagrangian
1909       call max_accel
1910       if (amax*d_time .gt. dvmax) then
1911         d_time=d_time*dvmax/amax
1912         if(me.eq.king.or..not.out1file) write (iout,*)
1913      &   "Time step reduced to",d_time,
1914      &   " because of too large initial acceleration."
1915       endif
1916       if(me.eq.king.or..not.out1file)then
1917        write(iout,*) "Potential energy and its components"
1918        call enerprint(potEcomp)
1919 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1920       endif
1921       potE=potEcomp(0)-potEcomp(20)
1922       totE=EK+potE
1923       itime=0
1924       if (ntwe.ne.0) call statout(itime)
1925       if(me.eq.king.or..not.out1file)
1926      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1927      &   " Kinetic energy",EK," potential energy",potE,
1928      &   " total energy",totE," maximum acceleration ",
1929      &   amax
1930       if (large) then
1931         write (iout,*) "Initial coordinates"
1932         do i=1,nres
1933           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1934      &    (c(j,i+nres),j=1,3)
1935         enddo
1936         write (iout,*) "Initial dC"
1937         do i=0,nres
1938           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1939      &    (dc(j,i+nres),j=1,3)
1940         enddo
1941         write (iout,*) "Initial velocities"
1942         write (iout,"(13x,' backbone ',23x,' side chain')")
1943         do i=0,nres
1944           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1945      &    (d_t(j,i+nres),j=1,3)
1946         enddo
1947         write (iout,*) "Initial accelerations"
1948         do i=0,nres
1949 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1950           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1951      &    (d_a(j,i+nres),j=1,3)
1952         enddo
1953       endif
1954       do i=0,2*nres
1955         do j=1,3
1956           dc_old(j,i)=dc(j,i)
1957           d_t_old(j,i)=d_t(j,i)
1958           d_a_old(j,i)=d_a(j,i)
1959         enddo
1960 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1961       enddo
1962       if (RESPA) then
1963 #ifdef MPI
1964       tt0 =MPI_Wtime()
1965 #else
1966       tt0 = tcpu()
1967 #endif
1968         call zerograd
1969         call etotal_short(energia_short)
1970 #ifdef TIMING_ENE
1971 #ifdef MPI
1972         t_eshort=t_eshort+MPI_Wtime()-tt0
1973 #else
1974         t_eshort=t_eshort+tcpu()-tt0
1975 #endif
1976 #endif
1977
1978         if(tnp .or. tnp1) then
1979          E_short=energia_short(0)
1980          HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
1981          Csplit=Hnose1
1982 c         Csplit =110
1983 c_new_var_csplit          Csplit=H0-E_long
1984 c          Csplit = H0-energia_short(0)
1985           write(iout,*) 'Csplit= ',Csplit
1986         endif
1987
1988
1989         call cartgrad
1990         call lagrangian
1991         if(.not.out1file .and. large) then
1992           write (iout,*) "energia_long",energia_long(0),
1993      &     " energia_short",energia_short(0),
1994      &     " total",energia_long(0)+energia_short(0)
1995           write (iout,*) "Initial fast-force accelerations"
1996           do i=0,nres
1997             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1998      &      (d_a(j,i+nres),j=1,3)
1999           enddo
2000         endif
2001 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2002         do i=0,2*nres
2003           do j=1,3
2004             d_a_short(j,i)=d_a(j,i)
2005           enddo
2006         enddo
2007 #ifdef MPI
2008         tt0=MPI_Wtime()
2009 #else
2010         tt0=tcpu()
2011 #endif
2012         call zerograd
2013         call etotal_long(energia_long)
2014 #ifdef TIMING_ENE
2015 #ifdef MPI
2016         t_elong=t_elong+MPI_Wtime()-tt0
2017 #else
2018         t_elong=t_elong+tcpu()-tt0
2019 #endif
2020 #endif
2021         call cartgrad
2022         call lagrangian
2023         if(.not.out1file .and. large) then
2024           write (iout,*) "energia_long",energia_long(0)
2025           write (iout,*) "Initial slow-force accelerations"
2026           do i=0,nres
2027             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2028      &      (d_a(j,i+nres),j=1,3)
2029           enddo
2030         endif
2031 #ifdef MPI
2032         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2033 #else
2034         t_enegrad=t_enegrad+tcpu()-tt0
2035 #endif
2036       endif
2037
2038
2039
2040       return
2041       end
2042 c-----------------------------------------------------------
2043       subroutine random_vel
2044       implicit real*8 (a-h,o-z)
2045       include 'DIMENSIONS'
2046       include 'COMMON.CONTROL'
2047       include 'COMMON.VAR'
2048       include 'COMMON.MD'
2049 #ifndef LANG0
2050       include 'COMMON.LANGEVIN'
2051 #else
2052       include 'COMMON.LANGEVIN.lang0'
2053 #endif
2054       include 'COMMON.CHAIN'
2055       include 'COMMON.DERIV'
2056       include 'COMMON.GEO'
2057       include 'COMMON.LOCAL'
2058       include 'COMMON.INTERACT'
2059       include 'COMMON.IOUNITS'
2060       include 'COMMON.NAMES'
2061       include 'COMMON.TIME1'
2062       double precision xv,sigv,lowb,highb
2063 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2064 c First generate velocities in the eigenspace of the G matrix
2065 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2066 c      call flush(iout)
2067 c      write (iout,*) "RANDOM_VEL dimen",dimen
2068       xv=0.0d0
2069       ii=0
2070       do i=1,dimen
2071         do k=1,3
2072           ii=ii+1
2073           sigv=dsqrt((Rb*t_bath)/geigen(i))
2074           lowb=-5*sigv
2075           highb=5*sigv
2076           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2077 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2078 c     &      " d_t_work_new",d_t_work_new(ii)
2079         enddo
2080       enddo
2081 c diagnostics
2082 c      Ek1=0.0d0
2083 c      ii=0
2084 c      do i=1,dimen
2085 c        do k=1,3
2086 c          ii=ii+1
2087 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2088 c        enddo
2089 c      enddo
2090 c      write (iout,*) "Ek from eigenvectors",Ek1
2091 c end diagnostics
2092 c Transform velocities to UNRES coordinate space
2093       do k=0,2
2094         do i=1,dimen
2095           ind=(i-1)*3+k+1
2096           d_t_work(ind)=0.0d0
2097           do j=1,dimen
2098             d_t_work(ind)=d_t_work(ind)
2099      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2100           enddo
2101 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2102 c          call flush(iout)
2103         enddo
2104       enddo
2105 c Transfer to the d_t vector
2106       do j=1,3
2107         d_t(j,0)=d_t_work(j)
2108       enddo
2109       ind=3
2110       do i=nnt,nct-1
2111         do j=1,3
2112           ind=ind+1
2113           d_t(j,i)=d_t_work(ind)
2114         enddo
2115       enddo
2116 c      do i=0,nres-1
2117 c        write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
2118 c      enddo
2119       do i=nnt,nct
2120         if (itype(i).ne.10) then
2121           do j=1,3
2122             ind=ind+1
2123             d_t(j,i+nres)=d_t_work(ind)
2124           enddo
2125         endif
2126       enddo
2127 c      call kinetic(EK)
2128 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2129 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2130 c      call flush(iout)
2131       return
2132       end
2133 #ifndef LANG0
2134 c-----------------------------------------------------------
2135       subroutine sd_verlet_p_setup
2136 c Sets up the parameters of stochastic Verlet algorithm
2137       implicit real*8 (a-h,o-z)
2138       include 'DIMENSIONS'
2139 #ifdef MPI
2140       include 'mpif.h'
2141 #endif
2142       include 'COMMON.CONTROL'
2143       include 'COMMON.VAR'
2144       include 'COMMON.MD'
2145 #ifndef LANG0
2146       include 'COMMON.LANGEVIN'
2147 #else
2148       include 'COMMON.LANGEVIN.lang0'
2149 #endif
2150       include 'COMMON.CHAIN'
2151       include 'COMMON.DERIV'
2152       include 'COMMON.GEO'
2153       include 'COMMON.LOCAL'
2154       include 'COMMON.INTERACT'
2155       include 'COMMON.IOUNITS'
2156       include 'COMMON.NAMES'
2157       include 'COMMON.TIME1'
2158       double precision emgdt(MAXRES6),
2159      & pterm,vterm,rho,rhoc,vsig,
2160      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2161      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2162      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2163       logical lprn /.false./
2164       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2165       double precision ktm
2166 #ifdef MPI
2167       tt0 = MPI_Wtime()
2168 #else
2169       tt0 = tcpu()
2170 #endif
2171 c
2172 c AL 8/17/04 Code adapted from tinker
2173 c
2174 c Get the frictional and random terms for stochastic dynamics in the
2175 c eigenspace of mass-scaled UNRES friction matrix
2176 c
2177       do i = 1, dimen
2178             gdt = fricgam(i) * d_time
2179 c
2180 c Stochastic dynamics reduces to simple MD for zero friction
2181 c
2182             if (gdt .le. zero) then
2183                pfric_vec(i) = 1.0d0
2184                vfric_vec(i) = d_time
2185                afric_vec(i) = 0.5d0 * d_time * d_time
2186                prand_vec(i) = 0.0d0
2187                vrand_vec1(i) = 0.0d0
2188                vrand_vec2(i) = 0.0d0
2189 c
2190 c Analytical expressions when friction coefficient is large
2191 c
2192             else
2193                if (gdt .ge. gdt_radius) then
2194                   egdt = dexp(-gdt)
2195                   pfric_vec(i) = egdt
2196                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2197                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2198                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2199                   vterm = 1.0d0 - egdt**2
2200                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2201 c
2202 c Use series expansions when friction coefficient is small
2203 c
2204                else
2205                   gdt2 = gdt * gdt
2206                   gdt3 = gdt * gdt2
2207                   gdt4 = gdt2 * gdt2
2208                   gdt5 = gdt2 * gdt3
2209                   gdt6 = gdt3 * gdt3
2210                   gdt7 = gdt3 * gdt4
2211                   gdt8 = gdt4 * gdt4
2212                   gdt9 = gdt4 * gdt5
2213                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2214      &                          - gdt5/120.0d0 + gdt6/720.0d0
2215      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2216      &                          - gdt9/362880.0d0) / fricgam(i)**2
2217                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2218                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2219                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2220      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2221      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2222      &                       + 127.0d0*gdt9/90720.0d0
2223                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2224      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2225      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2226      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2227                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2228      &                       - 17.0d0*gdt2/1280.0d0
2229      &                       + 17.0d0*gdt3/6144.0d0
2230      &                       + 40967.0d0*gdt4/34406400.0d0
2231      &                       - 57203.0d0*gdt5/275251200.0d0
2232      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2233                end if
2234 c
2235 c Compute the scaling factors of random terms for the nonzero friction case
2236 c
2237                ktm = 0.5d0*d_time/fricgam(i)
2238                psig = dsqrt(ktm*pterm) / fricgam(i)
2239                vsig = dsqrt(ktm*vterm)
2240                rhoc = dsqrt(1.0d0 - rho*rho)
2241                prand_vec(i) = psig
2242                vrand_vec1(i) = vsig * rho
2243                vrand_vec2(i) = vsig * rhoc
2244             end if
2245       end do
2246       if (lprn) then
2247       write (iout,*)
2248      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2249      &  " vrand_vec2"
2250       do i=1,dimen
2251         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2252      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2253       enddo
2254       endif
2255 c
2256 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2257 c
2258 #ifndef   LANG0
2259       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2260       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2261       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2262       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2263       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2264       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2265 #endif
2266 #ifdef MPI
2267       t_sdsetup=t_sdsetup+MPI_Wtime()
2268 #else
2269       t_sdsetup=t_sdsetup+tcpu()-tt0
2270 #endif
2271       return
2272       end
2273 c-------------------------------------------------------------
2274       subroutine eigtransf1(n,ndim,ab,d,c)
2275       implicit none
2276       integer n,ndim
2277       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2278       integer i,j,k
2279       do i=1,n
2280         do j=1,n
2281           c(i,j)=0.0d0
2282           do k=1,n
2283             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2284           enddo
2285         enddo
2286       enddo
2287       return
2288       end
2289 c-------------------------------------------------------------
2290       subroutine eigtransf(n,ndim,a,b,d,c)
2291       implicit none
2292       integer n,ndim
2293       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2294       integer i,j,k
2295       do i=1,n
2296         do j=1,n
2297           c(i,j)=0.0d0
2298           do k=1,n
2299             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2300           enddo
2301         enddo
2302       enddo
2303       return
2304       end
2305 c-------------------------------------------------------------
2306       subroutine sd_verlet1
2307 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2308       implicit real*8 (a-h,o-z)
2309       include 'DIMENSIONS'
2310       include 'COMMON.CONTROL'
2311       include 'COMMON.VAR'
2312       include 'COMMON.MD'
2313 #ifndef LANG0
2314       include 'COMMON.LANGEVIN'
2315 #else
2316       include 'COMMON.LANGEVIN.lang0'
2317 #endif
2318       include 'COMMON.CHAIN'
2319       include 'COMMON.DERIV'
2320       include 'COMMON.GEO'
2321       include 'COMMON.LOCAL'
2322       include 'COMMON.INTERACT'
2323       include 'COMMON.IOUNITS'
2324       include 'COMMON.NAMES'
2325       double precision stochforcvec(MAXRES6)
2326       common /stochcalc/ stochforcvec
2327       logical lprn /.false./
2328
2329 c      write (iout,*) "dc_old"
2330 c      do i=0,nres
2331 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2332 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2333 c      enddo
2334       do j=1,3
2335         dc_work(j)=dc_old(j,0)
2336         d_t_work(j)=d_t_old(j,0)
2337         d_a_work(j)=d_a_old(j,0)
2338       enddo
2339       ind=3
2340       do i=nnt,nct-1
2341         do j=1,3
2342           dc_work(ind+j)=dc_old(j,i)
2343           d_t_work(ind+j)=d_t_old(j,i)
2344           d_a_work(ind+j)=d_a_old(j,i)
2345         enddo
2346         ind=ind+3
2347       enddo
2348       do i=nnt,nct
2349         if (itype(i).ne.10) then
2350           do j=1,3
2351             dc_work(ind+j)=dc_old(j,i+nres)
2352             d_t_work(ind+j)=d_t_old(j,i+nres)
2353             d_a_work(ind+j)=d_a_old(j,i+nres)
2354           enddo
2355           ind=ind+3
2356         endif
2357       enddo
2358 #ifndef LANG0
2359       if (lprn) then
2360       write (iout,*)
2361      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2362      &  " vrand_mat2"
2363       do i=1,dimen
2364         do j=1,dimen
2365           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2366      &      vfric_mat(i,j),afric_mat(i,j),
2367      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2368         enddo
2369       enddo
2370       endif
2371       do i=1,dimen
2372         ddt1=0.0d0
2373         ddt2=0.0d0
2374         do j=1,dimen
2375           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2376      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2377           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2378           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2379         enddo
2380         d_t_work_new(i)=ddt1+0.5d0*ddt2
2381         d_t_work(i)=ddt1+ddt2
2382       enddo
2383 #endif
2384       do j=1,3
2385         dc(j,0)=dc_work(j)
2386         d_t(j,0)=d_t_work(j)
2387       enddo
2388       ind=3
2389       do i=nnt,nct-1
2390         do j=1,3
2391           dc(j,i)=dc_work(ind+j)
2392           d_t(j,i)=d_t_work(ind+j)
2393         enddo
2394         ind=ind+3
2395       enddo
2396       do i=nnt,nct
2397         if (itype(i).ne.10) then
2398           inres=i+nres
2399           do j=1,3
2400             dc(j,inres)=dc_work(ind+j)
2401             d_t(j,inres)=d_t_work(ind+j)
2402           enddo
2403           ind=ind+3
2404         endif
2405       enddo
2406       return
2407       end
2408 c--------------------------------------------------------------------------
2409       subroutine sd_verlet2
2410 c  Calculating the adjusted velocities for accelerations
2411       implicit real*8 (a-h,o-z)
2412       include 'DIMENSIONS'
2413       include 'COMMON.CONTROL'
2414       include 'COMMON.VAR'
2415       include 'COMMON.MD'
2416 #ifndef LANG0
2417       include 'COMMON.LANGEVIN'
2418 #else
2419       include 'COMMON.LANGEVIN.lang0'
2420 #endif
2421       include 'COMMON.CHAIN'
2422       include 'COMMON.DERIV'
2423       include 'COMMON.GEO'
2424       include 'COMMON.LOCAL'
2425       include 'COMMON.INTERACT'
2426       include 'COMMON.IOUNITS'
2427       include 'COMMON.NAMES'
2428       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2429       common /stochcalc/ stochforcvec
2430 c
2431 c Compute the stochastic forces which contribute to velocity change
2432 c
2433       call stochastic_force(stochforcvecV)
2434
2435 #ifndef LANG0
2436       do i=1,dimen
2437         ddt1=0.0d0
2438         ddt2=0.0d0
2439         do j=1,dimen
2440           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2441           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2442      &     vrand_mat2(i,j)*stochforcvecV(j)
2443         enddo
2444         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2445       enddo
2446 #endif
2447       do j=1,3
2448         d_t(j,0)=d_t_work(j)
2449       enddo
2450       ind=3
2451       do i=nnt,nct-1
2452         do j=1,3
2453           d_t(j,i)=d_t_work(ind+j)
2454         enddo
2455         ind=ind+3
2456       enddo
2457       do i=nnt,nct
2458         if (itype(i).ne.10) then
2459           inres=i+nres
2460           do j=1,3
2461             d_t(j,inres)=d_t_work(ind+j)
2462           enddo
2463           ind=ind+3
2464         endif
2465       enddo
2466       return
2467       end
2468 c-----------------------------------------------------------
2469       subroutine sd_verlet_ciccotti_setup
2470 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2471 c version
2472       implicit real*8 (a-h,o-z)
2473       include 'DIMENSIONS'
2474 #ifdef MPI
2475       include 'mpif.h'
2476 #endif
2477       include 'COMMON.CONTROL'
2478       include 'COMMON.VAR'
2479       include 'COMMON.MD'
2480 #ifndef LANG0
2481       include 'COMMON.LANGEVIN'
2482 #else
2483       include 'COMMON.LANGEVIN.lang0'
2484 #endif
2485       include 'COMMON.CHAIN'
2486       include 'COMMON.DERIV'
2487       include 'COMMON.GEO'
2488       include 'COMMON.LOCAL'
2489       include 'COMMON.INTERACT'
2490       include 'COMMON.IOUNITS'
2491       include 'COMMON.NAMES'
2492       include 'COMMON.TIME1'
2493       double precision emgdt(MAXRES6),
2494      & pterm,vterm,rho,rhoc,vsig,
2495      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2496      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2497      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2498       logical lprn /.false./
2499       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2500       double precision ktm
2501 #ifdef MPI
2502       tt0 = MPI_Wtime()
2503 #else
2504       tt0 = tcpu()
2505 #endif
2506 c
2507 c AL 8/17/04 Code adapted from tinker
2508 c
2509 c Get the frictional and random terms for stochastic dynamics in the
2510 c eigenspace of mass-scaled UNRES friction matrix
2511 c
2512       do i = 1, dimen
2513             write (iout,*) "i",i," fricgam",fricgam(i)
2514             gdt = fricgam(i) * d_time
2515 c
2516 c Stochastic dynamics reduces to simple MD for zero friction
2517 c
2518             if (gdt .le. zero) then
2519                pfric_vec(i) = 1.0d0
2520                vfric_vec(i) = d_time
2521                afric_vec(i) = 0.5d0*d_time*d_time
2522                prand_vec(i) = afric_vec(i)
2523                vrand_vec2(i) = vfric_vec(i)
2524 c
2525 c Analytical expressions when friction coefficient is large
2526 c
2527             else
2528                egdt = dexp(-gdt)
2529                pfric_vec(i) = egdt
2530                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2531                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2532                prand_vec(i) = afric_vec(i)
2533                vrand_vec2(i) = vfric_vec(i)
2534 c
2535 c Compute the scaling factors of random terms for the nonzero friction case
2536 c
2537 c               ktm = 0.5d0*d_time/fricgam(i)
2538 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2539 c               vsig = dsqrt(ktm*vterm)
2540 c               prand_vec(i) = psig*afric_vec(i)
2541 c               vrand_vec2(i) = vsig*vfric_vec(i)
2542             end if
2543       end do
2544       if (lprn) then
2545       write (iout,*)
2546      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2547      &  " vrand_vec2"
2548       do i=1,dimen
2549         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2550      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2551       enddo
2552       endif
2553 c
2554 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2555 c
2556       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2557       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2558       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2559       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2560       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2561 #ifdef MPI
2562       t_sdsetup=t_sdsetup+MPI_Wtime()
2563 #else
2564       t_sdsetup=t_sdsetup+tcpu()-tt0
2565 #endif
2566       return
2567       end
2568 c-------------------------------------------------------------
2569       subroutine sd_verlet1_ciccotti
2570 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2571       implicit real*8 (a-h,o-z)
2572       include 'DIMENSIONS'
2573 #ifdef MPI
2574       include 'mpif.h'
2575 #endif
2576       include 'COMMON.CONTROL'
2577       include 'COMMON.VAR'
2578       include 'COMMON.MD'
2579 #ifndef LANG0
2580       include 'COMMON.LANGEVIN'
2581 #else
2582       include 'COMMON.LANGEVIN.lang0'
2583 #endif
2584       include 'COMMON.CHAIN'
2585       include 'COMMON.DERIV'
2586       include 'COMMON.GEO'
2587       include 'COMMON.LOCAL'
2588       include 'COMMON.INTERACT'
2589       include 'COMMON.IOUNITS'
2590       include 'COMMON.NAMES'
2591       double precision stochforcvec(MAXRES6)
2592       common /stochcalc/ stochforcvec
2593       logical lprn /.false./
2594
2595 c      write (iout,*) "dc_old"
2596 c      do i=0,nres
2597 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2598 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2599 c      enddo
2600       do j=1,3
2601         dc_work(j)=dc_old(j,0)
2602         d_t_work(j)=d_t_old(j,0)
2603         d_a_work(j)=d_a_old(j,0)
2604       enddo
2605       ind=3
2606       do i=nnt,nct-1
2607         do j=1,3
2608           dc_work(ind+j)=dc_old(j,i)
2609           d_t_work(ind+j)=d_t_old(j,i)
2610           d_a_work(ind+j)=d_a_old(j,i)
2611         enddo
2612         ind=ind+3
2613       enddo
2614       do i=nnt,nct
2615         if (itype(i).ne.10) then
2616           do j=1,3
2617             dc_work(ind+j)=dc_old(j,i+nres)
2618             d_t_work(ind+j)=d_t_old(j,i+nres)
2619             d_a_work(ind+j)=d_a_old(j,i+nres)
2620           enddo
2621           ind=ind+3
2622         endif
2623       enddo
2624
2625 #ifndef LANG0
2626       if (lprn) then
2627       write (iout,*)
2628      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2629      &  " vrand_mat2"
2630       do i=1,dimen
2631         do j=1,dimen
2632           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2633      &      vfric_mat(i,j),afric_mat(i,j),
2634      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2635         enddo
2636       enddo
2637       endif
2638       do i=1,dimen
2639         ddt1=0.0d0
2640         ddt2=0.0d0
2641         do j=1,dimen
2642           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2643      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2644           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2645           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2646         enddo
2647         d_t_work_new(i)=ddt1+0.5d0*ddt2
2648         d_t_work(i)=ddt1+ddt2
2649       enddo
2650 #endif
2651       do j=1,3
2652         dc(j,0)=dc_work(j)
2653         d_t(j,0)=d_t_work(j)
2654       enddo
2655       ind=3
2656       do i=nnt,nct-1
2657         do j=1,3
2658           dc(j,i)=dc_work(ind+j)
2659           d_t(j,i)=d_t_work(ind+j)
2660         enddo
2661         ind=ind+3
2662       enddo
2663       do i=nnt,nct
2664         if (itype(i).ne.10) then
2665           inres=i+nres
2666           do j=1,3
2667             dc(j,inres)=dc_work(ind+j)
2668             d_t(j,inres)=d_t_work(ind+j)
2669           enddo
2670           ind=ind+3
2671         endif
2672       enddo
2673       return
2674       end
2675 c--------------------------------------------------------------------------
2676       subroutine sd_verlet2_ciccotti
2677 c  Calculating the adjusted velocities for accelerations
2678       implicit real*8 (a-h,o-z)
2679       include 'DIMENSIONS'
2680       include 'COMMON.CONTROL'
2681       include 'COMMON.VAR'
2682       include 'COMMON.MD'
2683 #ifndef LANG0
2684       include 'COMMON.LANGEVIN'
2685 #else
2686       include 'COMMON.LANGEVIN.lang0'
2687 #endif
2688       include 'COMMON.CHAIN'
2689       include 'COMMON.DERIV'
2690       include 'COMMON.GEO'
2691       include 'COMMON.LOCAL'
2692       include 'COMMON.INTERACT'
2693       include 'COMMON.IOUNITS'
2694       include 'COMMON.NAMES'
2695       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2696       common /stochcalc/ stochforcvec
2697 c
2698 c Compute the stochastic forces which contribute to velocity change
2699 c
2700       call stochastic_force(stochforcvecV)
2701 #ifndef LANG0
2702       do i=1,dimen
2703         ddt1=0.0d0
2704         ddt2=0.0d0
2705         do j=1,dimen
2706
2707           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2708 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2709           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2710         enddo
2711         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2712       enddo
2713 #endif
2714       do j=1,3
2715         d_t(j,0)=d_t_work(j)
2716       enddo
2717       ind=3
2718       do i=nnt,nct-1
2719         do j=1,3
2720           d_t(j,i)=d_t_work(ind+j)
2721         enddo
2722         ind=ind+3
2723       enddo
2724       do i=nnt,nct
2725         if (itype(i).ne.10) then
2726           inres=i+nres
2727           do j=1,3
2728             d_t(j,inres)=d_t_work(ind+j)
2729           enddo
2730           ind=ind+3
2731         endif
2732       enddo
2733       return
2734       end
2735 #endif
2736 c------------------------------------------------------
2737       double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
2738       implicit none
2739       double precision ek,s,e,pi,Q,t_bath,Rb
2740       integer dimenl
2741       Rb=0.001986d0
2742       HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
2743 c      print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
2744 c     &      pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
2745       return
2746       end
2747 c-----------------------------------------------------------------
2748       double precision function HNose_nh(eki,e)
2749       implicit real*8 (a-h,o-z)
2750       include 'DIMENSIONS'
2751       include 'COMMON.MD'
2752       HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
2753       do i=2,nnos
2754         HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
2755       enddo
2756 c      write(4,'(5e15.5)')
2757 c     &       vlogs(1),xlogs(1),HNose,eki,e
2758       return
2759       end
2760 c-----------------------------------------------------------------
2761       SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
2762       implicit real*8 (a-h,o-z)
2763       include 'DIMENSIONS'
2764       include 'COMMON.MD'
2765       double precision akin,gnkt,dt,aa,gkt,scale
2766       double precision wdti(maxyosh),wdti2(maxyosh),
2767      &                 wdti4(maxyosh),wdti8(maxyosh)
2768       integer i,iresn,iyosh,inos,nnos1
2769
2770       dt=d_time
2771       nnos1=nnos+1
2772       GKT = Rb*t_bath
2773       GNKT = dimen3*GKT
2774       akin=akin*2
2775
2776
2777 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
2778 C INTEGRATION FROM t=0 TO t=DT/2
2779 C GET THE TOTAL KINETIC ENERGY
2780       SCALE = 1.D0
2781 c      CALL GETKINP(MASS,VX,VY,VZ,AKIN)
2782 C UPDATE THE FORCES
2783       GLOGS(1) = (AKIN - GNKT)/QMASS(1)
2784 C START THE MULTIPLE TIME STEP PROCEDURE
2785       DO IRESN = 1,NRESN
2786        DO IYOSH = 1,NYOSH
2787 C UPDATE THE THERMOSTAT VELOCITIES
2788         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2789         DO INOS = 1,NNOS-1
2790          AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
2791          VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
2792      &          + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
2793         ENDDO
2794 C UPDATE THE PARTICLE VELOCITIES
2795         AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
2796         SCALE = SCALE*AA
2797 C UPDATE THE FORCES
2798         GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
2799 C UPDATE THE THERMOSTAT POSITIONS
2800         DO INOS = 1,NNOS
2801          XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
2802         ENDDO
2803 C UPDATE THE THERMOSTAT VELOCITIES
2804         DO INOS = 1,NNOS-1
2805          AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
2806          VLOGS(INOS) = VLOGS(INOS)*AA*AA
2807      &      + WDTI4(IYOSH)*GLOGS(INOS)*AA
2808          GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
2809      &      -GKT)/QMASS(INOS+1)
2810         ENDDO
2811         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
2812        ENDDO
2813       ENDDO
2814 C UPDATE THE PARTICLE VELOCITIES
2815 c outside of this subroutine
2816 c      DO I = 1,N
2817 c       VX(I) = VX(I)*SCALE
2818 c       VY(I) = VY(I)*SCALE
2819 c       VZ(I) = VZ(I)*SCALE
2820 c      ENDDO
2821       RETURN
2822       END
2823 c-----------------------------------------------------------------
2824       subroutine tnp1_respa_i_step1
2825 c Applying Nose-Poincare algorithm - step 1 to coordinates
2826 c JPSJ 70 75 (2001) S. Nose
2827 c
2828 c d_t is not updated here
2829 c
2830       implicit real*8 (a-h,o-z)
2831       include 'DIMENSIONS'
2832       include 'COMMON.CONTROL'
2833       include 'COMMON.VAR'
2834       include 'COMMON.MD'
2835       include 'COMMON.CHAIN'
2836       include 'COMMON.DERIV'
2837       include 'COMMON.GEO'
2838       include 'COMMON.LOCAL'
2839       include 'COMMON.INTERACT'
2840       include 'COMMON.IOUNITS'
2841       include 'COMMON.NAMES'
2842       double precision adt,adt2,tmp
2843
2844       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2845       s12_np=s_np*tmp**2
2846       pistar=pi_np/tmp
2847       s12_dt=d_time/s12_np
2848       d_time_s12=d_time*0.5*s12_np
2849
2850       do j=1,3
2851         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2852         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2853       enddo
2854       do i=nnt,nct-1
2855         do j=1,3
2856           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2857           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2858         enddo
2859       enddo
2860       do i=nnt,nct
2861         if (itype(i).ne.10) then
2862           inres=i+nres
2863           do j=1,3
2864            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2865            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2866           enddo
2867         endif
2868       enddo
2869       return
2870       end
2871 c---------------------------------------------------------------------
2872       subroutine tnp1_respa_i_step2
2873 c  Step 2 of the velocity Verlet algorithm: update velocities
2874       implicit real*8 (a-h,o-z)
2875       include 'DIMENSIONS'
2876       include 'COMMON.CONTROL'
2877       include 'COMMON.VAR'
2878       include 'COMMON.MD'
2879       include 'COMMON.CHAIN'
2880       include 'COMMON.DERIV'
2881       include 'COMMON.GEO'
2882       include 'COMMON.LOCAL'
2883       include 'COMMON.INTERACT'
2884       include 'COMMON.IOUNITS'
2885       include 'COMMON.NAMES'
2886
2887       double precision d_time_s12
2888
2889       do i=0,2*nres
2890        do j=1,3
2891         d_t(j,i)=d_t_new(j,i)
2892        enddo
2893       enddo
2894
2895       call kinetic(EK)
2896       EK=EK/s12_np**2
2897
2898       d_time_s12=0.5d0*s12_np*d_time
2899
2900       do j=1,3
2901         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
2902       enddo
2903       do i=nnt,nct-1
2904         do j=1,3
2905           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
2906         enddo
2907       enddo
2908       do i=nnt,nct
2909         if (itype(i).ne.10) then
2910           inres=i+nres
2911           do j=1,3
2912             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
2913           enddo
2914         endif
2915       enddo
2916
2917       pistar=pistar+(EK-0.5*(E_old+potE)
2918      &       -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
2919       tmp=1+pistar/(2*Q_np)*0.5*d_time
2920       s_np=s12_np*tmp**2
2921       pi_np=pistar/tmp
2922
2923       return
2924       end
2925 c-------------------------------------------------------
2926
2927       subroutine tnp1_step1
2928 c Applying Nose-Poincare algorithm - step 1 to coordinates
2929 c JPSJ 70 75 (2001) S. Nose
2930 c
2931 c d_t is not updated here
2932 c
2933       implicit real*8 (a-h,o-z)
2934       include 'DIMENSIONS'
2935       include 'COMMON.CONTROL'
2936       include 'COMMON.VAR'
2937       include 'COMMON.MD'
2938       include 'COMMON.CHAIN'
2939       include 'COMMON.DERIV'
2940       include 'COMMON.GEO'
2941       include 'COMMON.LOCAL'
2942       include 'COMMON.INTERACT'
2943       include 'COMMON.IOUNITS'
2944       include 'COMMON.NAMES'
2945       double precision adt,adt2,tmp
2946
2947       tmp=1+pi_np/(2*Q_np)*0.5*d_time
2948       s12_np=s_np*tmp**2
2949       pistar=pi_np/tmp
2950       s12_dt=d_time/s12_np
2951       d_time_s12=d_time*0.5*s12_np
2952
2953       do j=1,3
2954         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
2955         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
2956       enddo
2957       do i=nnt,nct-1
2958         do j=1,3
2959           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
2960           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
2961         enddo
2962       enddo
2963       do i=nnt,nct
2964         if (itype(i).ne.10) then
2965           inres=i+nres
2966           do j=1,3
2967            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
2968            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
2969           enddo
2970         endif
2971       enddo
2972       return
2973       end
2974 c---------------------------------------------------------------------
2975       subroutine tnp1_step2
2976 c  Step 2 of the velocity Verlet algorithm: update velocities
2977       implicit real*8 (a-h,o-z)
2978       include 'DIMENSIONS'
2979       include 'COMMON.CONTROL'
2980       include 'COMMON.VAR'
2981       include 'COMMON.MD'
2982       include 'COMMON.CHAIN'
2983       include 'COMMON.DERIV'
2984       include 'COMMON.GEO'
2985       include 'COMMON.LOCAL'
2986       include 'COMMON.INTERACT'
2987       include 'COMMON.IOUNITS'
2988       include 'COMMON.NAMES'
2989
2990       double precision d_time_s12
2991
2992       do i=0,2*nres
2993        do j=1,3
2994         d_t(j,i)=d_t_new(j,i)
2995        enddo
2996       enddo
2997
2998       call kinetic(EK)
2999       EK=EK/s12_np**2
3000
3001       d_time_s12=0.5d0*s12_np*d_time
3002
3003       do j=1,3
3004         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
3005       enddo
3006       do i=nnt,nct-1
3007         do j=1,3
3008           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
3009         enddo
3010       enddo
3011       do i=nnt,nct
3012         if (itype(i).ne.10) then
3013           inres=i+nres
3014           do j=1,3
3015             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
3016           enddo
3017         endif
3018       enddo
3019
3020 cd      write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
3021       pistar=pistar+(EK-0.5*(E_old+potE)
3022      &       -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
3023       tmp=1+pistar/(2*Q_np)*0.5*d_time
3024       s_np=s12_np*tmp**2
3025       pi_np=pistar/tmp
3026
3027       return
3028       end
3029
3030 c-----------------------------------------------------------------
3031       subroutine tnp_respa_i_step1
3032 c Applying Nose-Poincare algorithm - step 1 to coordinates
3033 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3034 c
3035 c d_t is not updated here, it is destroyed
3036 c
3037       implicit real*8 (a-h,o-z)
3038       include 'DIMENSIONS'
3039       include 'COMMON.CONTROL'
3040       include 'COMMON.VAR'
3041       include 'COMMON.MD'
3042       include 'COMMON.CHAIN'
3043       include 'COMMON.DERIV'
3044       include 'COMMON.GEO'
3045       include 'COMMON.LOCAL'
3046       include 'COMMON.INTERACT'
3047       include 'COMMON.IOUNITS'
3048       include 'COMMON.NAMES'
3049       double precision C_np,d_time_s,tmp,d_time_ss
3050
3051       d_time_s=d_time*0.5*s_np
3052 ct2      d_time_s=d_time*0.5*s12_np
3053
3054       do j=1,3
3055         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3056       enddo
3057       do i=nnt,nct-1
3058         do j=1,3
3059           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3060         enddo
3061       enddo
3062       do i=nnt,nct
3063         if (itype(i).ne.10) then
3064           inres=i+nres
3065           do j=1,3
3066            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3067           enddo
3068         endif
3069       enddo
3070
3071       do i=0,2*nres
3072        do j=1,3
3073         d_t(j,i)=d_t_new(j,i)
3074        enddo
3075       enddo
3076
3077       call kinetic(EK)
3078       EK=EK/s_np**2
3079
3080       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
3081      &                     -pi_np
3082
3083       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3084       tmp=0.5*d_time*pistar/Q_np
3085       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3086
3087       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3088 ct2      d_time_ss=d_time/s12_np
3089 c      d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
3090
3091       do j=1,3
3092         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3093       enddo
3094       do i=nnt,nct-1
3095         do j=1,3
3096           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3097         enddo
3098       enddo
3099       do i=nnt,nct
3100         if (itype(i).ne.10) then
3101           inres=i+nres
3102           do j=1,3
3103            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3104           enddo
3105         endif
3106       enddo
3107
3108       return
3109       end
3110 c---------------------------------------------------------------------
3111
3112       subroutine tnp_respa_i_step2
3113 c  Step 2 of the velocity Verlet algorithm: update velocities
3114       implicit real*8 (a-h,o-z)
3115       include 'DIMENSIONS'
3116       include 'COMMON.CONTROL'
3117       include 'COMMON.VAR'
3118       include 'COMMON.MD'
3119       include 'COMMON.CHAIN'
3120       include 'COMMON.DERIV'
3121       include 'COMMON.GEO'
3122       include 'COMMON.LOCAL'
3123       include 'COMMON.INTERACT'
3124       include 'COMMON.IOUNITS'
3125       include 'COMMON.NAMES'
3126
3127       double precision d_time_s
3128
3129       EK=EK*(s_np/s12_np)**2
3130       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3131       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
3132      &                              -HNose1+Csplit)
3133
3134 cr      print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
3135       d_time_s=d_time*0.5*s12_np
3136 c      d_time_s=d_time*0.5*s_np
3137
3138       do j=1,3
3139         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3140       enddo
3141       do i=nnt,nct-1
3142         do j=1,3
3143           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3144         enddo
3145       enddo
3146       do i=nnt,nct
3147         if (itype(i).ne.10) then
3148           inres=i+nres
3149           do j=1,3
3150             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3151           enddo
3152         endif
3153       enddo
3154
3155       s_np=s12_np
3156
3157       return
3158       end
3159 c-----------------------------------------------------------------
3160       subroutine tnp_respa_step1
3161 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
3162 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3163 c
3164 c d_t is not updated here, it is destroyed
3165 c
3166       implicit real*8 (a-h,o-z)
3167       include 'DIMENSIONS'
3168       include 'COMMON.CONTROL'
3169       include 'COMMON.VAR'
3170       include 'COMMON.MD'
3171       include 'COMMON.CHAIN'
3172       include 'COMMON.DERIV'
3173       include 'COMMON.GEO'
3174       include 'COMMON.LOCAL'
3175       include 'COMMON.INTERACT'
3176       include 'COMMON.IOUNITS'
3177       include 'COMMON.NAMES'
3178       double precision C_np,d_time_s,tmp,d_time_ss
3179       double precision energia(0:n_ene)
3180
3181       d_time_s=d_time*0.5*s_np
3182
3183       do j=1,3
3184         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3185       enddo
3186       do i=nnt,nct-1
3187         do j=1,3
3188           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3189         enddo
3190       enddo
3191       do i=nnt,nct
3192         if (itype(i).ne.10) then
3193           inres=i+nres
3194           do j=1,3
3195            d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3196           enddo
3197         endif
3198       enddo
3199
3200
3201 c      C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3202 c     &                     -pi_np
3203 c
3204 c      pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3205 c      tmp=0.5*d_time*pistar/Q_np
3206 c      s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3207 c      write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3208
3209 ct1      pi_np=pistar
3210 c      sold_np=s_np
3211 c      s_np=s12_np
3212
3213 c-------------------------------------
3214 c test of reviewer's comment
3215        pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3216 cr       print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
3217 c-------------------------------------
3218
3219       return
3220       end
3221 c---------------------------------------------------------------------
3222       subroutine tnp_respa_step2
3223 c  Step 2 of the velocity Verlet algorithm: update velocities for RESPA
3224       implicit real*8 (a-h,o-z)
3225       include 'DIMENSIONS'
3226       include 'COMMON.CONTROL'
3227       include 'COMMON.VAR'
3228       include 'COMMON.MD'
3229       include 'COMMON.CHAIN'
3230       include 'COMMON.DERIV'
3231       include 'COMMON.GEO'
3232       include 'COMMON.LOCAL'
3233       include 'COMMON.INTERACT'
3234       include 'COMMON.IOUNITS'
3235       include 'COMMON.NAMES'
3236
3237       double precision d_time_s
3238
3239 ct1      s12_np=s_np
3240 ct2      pistar=pi_np
3241
3242 ct      call kinetic(EK)
3243 ct      HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3244 ct      pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3245 ct     &                              -0.5*d_time*(HNose1-H0)
3246
3247 c-------------------------------------
3248 c test of reviewer's comment
3249       pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
3250 cr      print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
3251 c-------------------------------------
3252       d_time_s=d_time*0.5*s_np
3253
3254       do j=1,3
3255         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
3256       enddo
3257       do i=nnt,nct-1
3258         do j=1,3
3259           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
3260         enddo
3261       enddo
3262       do i=nnt,nct
3263         if (itype(i).ne.10) then
3264           inres=i+nres
3265           do j=1,3
3266             d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
3267           enddo
3268         endif
3269       enddo
3270
3271 cd      s_np=s12_np
3272
3273       return
3274       end
3275 c---------------------------------------------------------------------
3276       subroutine tnp_step1
3277 c Applying Nose-Poincare algorithm - step 1 to coordinates
3278 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
3279 c
3280 c d_t is not updated here, it is destroyed
3281 c
3282       implicit real*8 (a-h,o-z)
3283       include 'DIMENSIONS'
3284       include 'COMMON.CONTROL'
3285       include 'COMMON.VAR'
3286       include 'COMMON.MD'
3287       include 'COMMON.CHAIN'
3288       include 'COMMON.DERIV'
3289       include 'COMMON.GEO'
3290       include 'COMMON.LOCAL'
3291       include 'COMMON.INTERACT'
3292       include 'COMMON.IOUNITS'
3293       include 'COMMON.NAMES'
3294       double precision C_np,d_time_s,tmp,d_time_ss
3295
3296       d_time_s=d_time*0.5*s_np
3297
3298       do j=1,3
3299         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
3300       enddo
3301       do i=nnt,nct-1
3302         do j=1,3
3303           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
3304         enddo
3305       enddo
3306       do i=nnt,nct
3307         if (itype(i).ne.10) then
3308           inres=i+nres
3309           do j=1,3
3310            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
3311           enddo
3312         endif
3313       enddo
3314
3315       do i=0,2*nres
3316        do j=1,3
3317         d_t(j,i)=d_t_new(j,i)
3318        enddo
3319       enddo
3320
3321       call kinetic(EK)
3322       EK=EK/s_np**2
3323
3324       C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
3325      &                     -pi_np
3326
3327       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
3328       tmp=0.5*d_time*pistar/Q_np
3329       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
3330 c      write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
3331
3332       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
3333
3334       do j=1,3
3335         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
3336       enddo
3337       do i=nnt,nct-1
3338         do j=1,3
3339           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
3340         enddo
3341       enddo
3342       do i=nnt,nct
3343         if (itype(i).ne.10) then
3344           inres=i+nres
3345           do j=1,3
3346            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
3347           enddo
3348         endif
3349       enddo
3350
3351       return
3352       end
3353 c-----------------------------------------------------------------
3354       subroutine tnp_step2
3355 c  Step 2 of the velocity Verlet algorithm: update velocities
3356       implicit real*8 (a-h,o-z)
3357       include 'DIMENSIONS'
3358       include 'COMMON.CONTROL'
3359       include 'COMMON.VAR'
3360       include 'COMMON.MD'
3361       include 'COMMON.CHAIN'
3362       include 'COMMON.DERIV'
3363       include 'COMMON.GEO'
3364       include 'COMMON.LOCAL'
3365       include 'COMMON.INTERACT'
3366       include 'COMMON.IOUNITS'
3367       include 'COMMON.NAMES'
3368
3369       double precision d_time_s
3370
3371       EK=EK*(s_np/s12_np)**2
3372       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
3373       pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
3374      &                              -0.5*d_time*(HNose1-H0)
3375
3376 cd      write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
3377       d_time_s=d_time*0.5*s12_np
3378
3379       do j=1,3
3380         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
3381       enddo
3382       do i=nnt,nct-1
3383         do j=1,3
3384           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
3385         enddo
3386       enddo
3387       do i=nnt,nct
3388         if (itype(i).ne.10) then
3389           inres=i+nres
3390           do j=1,3
3391             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
3392           enddo
3393         endif
3394       enddo
3395
3396       s_np=s12_np
3397
3398       return
3399       end
3400
3401       subroutine hmc_test(itime)
3402       implicit real*8 (a-h,o-z)
3403       include 'DIMENSIONS'
3404       include 'COMMON.CONTROL'
3405       include 'COMMON.MD'
3406       include 'COMMON.CHAIN'
3407
3408            hmc_acc=hmc_acc+1
3409            delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
3410            if (delta .lt. -50.0d0) then
3411                 delta=0.0d0
3412            else
3413                 delta=dexp(delta)
3414            endif
3415            xxx=ran_number(0.0d0,1.0d0)
3416
3417            if (me.eq.king .or. .not. out1file)
3418      &       write(iout,'(a8,i5,6f10.4)')
3419      &        'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
3420
3421            if (delta .le. xxx) then
3422             do i=1,2*nres
3423              do j=1,3
3424               dc(j,i)=dc_hmc(j,i)
3425              enddo
3426             enddo
3427             itime=itime-hmc
3428             totT=totThmc
3429            else
3430             if (me.eq.king .or. .not. out1file)
3431      &       write(iout,*) 'HMC accepting new'
3432             totThmc=totT
3433             do i=1,2*nres
3434              do j=1,3
3435               dc_hmc(j,i)=dc(j,i)
3436              enddo
3437             enddo
3438            endif
3439
3440            call chainbuild_cart
3441            call random_vel
3442            do i=0,2*nres
3443             do j=1,3
3444               d_t_old(j,i)=d_t(j,i)
3445             enddo
3446            enddo
3447            call kinetic(EK)
3448            kinetic_T=2.0d0/(dimen3*Rb)*EK
3449            call etotal(potEcomp)
3450            potE=potEcomp(0)
3451            hmc_etot=potE+EK
3452            if (me.eq.king .or. .not. out1file)
3453      &      write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK
3454
3455
3456       return
3457       end