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