source/unres/src_MD/MD_A-MTS.F

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