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