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