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