source/unres/src-5hdiag-tmp/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       include 'COMMON.LAGRANGE'
  16 #ifndef LANG0
  17       include 'COMMON.LANGEVIN'
  18 #else
  19       include 'COMMON.LANGEVIN.lang0'
  20 #endif
  21       include 'COMMON.CHAIN'
  22       include 'COMMON.DERIV'
  23       include 'COMMON.GEO'
  24       include 'COMMON.LOCAL'
  25       include 'COMMON.INTERACT'
  26       include 'COMMON.IOUNITS'
  27       include 'COMMON.NAMES'
  28       include 'COMMON.TIME1'
  29       include 'COMMON.HAIRPIN'
  30       double precision cm(3),L(3),vcm(3)
  31 #ifdef VOUT
  32       double precision v_work(maxres6),v_transf(maxres6)
  33 #endif
  34       integer ilen,rstcount
  35       external ilen
  36       character*50 tytul
  37       common /gucio/ cm
  38       integer itime
  39       logical ovrtim
  40 c
  41 #ifdef MPI
  42       if (ilen(tmpdir).gt.0)
  43      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
  44      &        //liczba(:ilen(liczba))//'.rst')
  45 #else
  46       if (ilen(tmpdir).gt.0)
  47      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
  48 #endif
  49       t_MDsetup=0.0d0
  50       t_langsetup=0.0d0
  51       t_MD=0.0d0
  52       t_enegrad=0.0d0
  53       t_sdsetup=0.0d0
  54       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  55 #ifdef MPI
  56       tt0=MPI_Wtime()
  57 #else
  58       tt0 = tcpu()
  59 #endif
  60 c Determine the inverse of the inertia matrix.
  61       call setup_MD_matrices
  62 c Initialize MD
  63       call init_MD
  64 #ifdef MPI
  65       t_MDsetup = MPI_Wtime()-tt0
  66 #else
  67       t_MDsetup = tcpu()-tt0
  68 #endif
  69       rstcount=0
  70 c   Entering the MD loop
  71 #ifdef MPI
  72       tt0 = MPI_Wtime()
  73 #else
  74       tt0 = tcpu()
  75 #endif
  76       if (lang.eq.2 .or. lang.eq.3) then
  77 #ifndef   LANG0
  78         call setup_fricmat
  79         if (lang.eq.2) then
  80           call sd_verlet_p_setup
  81         else
  82           call sd_verlet_ciccotti_setup
  83         endif
  84         do i=1,dimen
  85           do j=1,dimen
  86             pfric0_mat(i,j,0)=pfric_mat(i,j)
  87             afric0_mat(i,j,0)=afric_mat(i,j)
  88             vfric0_mat(i,j,0)=vfric_mat(i,j)
  89             prand0_mat(i,j,0)=prand_mat(i,j)
  90             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
  91             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
  92           enddo
  93         enddo
  94         flag_stoch(0)=.true.
  95         do i=1,maxflag_stoch
  96           flag_stoch(i)=.false.
  97         enddo
  98 #else
  99         write (iout,*)
 100      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 101 #ifdef MPI
 102         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 103 #endif
 104         stop
 105 #endif
 106       else if (lang.eq.1 .or. lang.eq.4) then
 107         call setup_fricmat
 108       endif
 109 #ifdef MPI
 110       t_langsetup=MPI_Wtime()-tt0
 111       tt0=MPI_Wtime()
 112 #else
 113       t_langsetup=tcpu()-tt0
 114       tt0=tcpu()
 115 #endif
 116       do itime=1,n_timestep
 117         if (ovrtim()) exit
 118         if (large.and. mod(itime,ntwe).eq.0)
 119      &    write (iout,*) "itime",itime
 120         rstcount=rstcount+1
 121         if (lang.gt.0 .and. surfarea .and.
 122      &      mod(itime,reset_fricmat).eq.0) then
 123           if (lang.eq.2 .or. lang.eq.3) then
 124 #ifndef LANG0
 125             call setup_fricmat
 126             if (lang.eq.2) then
 127               call sd_verlet_p_setup
 128             else
 129               call sd_verlet_ciccotti_setup
 130             endif
 131             do i=1,dimen
 132               do j=1,dimen
 133                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 134                 afric0_mat(i,j,0)=afric_mat(i,j)
 135                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 136                 prand0_mat(i,j,0)=prand_mat(i,j)
 137                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 138                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 139               enddo
 140             enddo
 141             flag_stoch(0)=.true.
 142             do i=1,maxflag_stoch
 143               flag_stoch(i)=.false.
 144             enddo
 145 #endif
 146           else if (lang.eq.1 .or. lang.eq.4) then
 147             call setup_fricmat
 148           endif
 149           write (iout,'(a,i10)')
 150      &      "Friction matrix reset based on surface area, itime",itime
 151         endif
 152         if (reset_vel .and. tbf .and. lang.eq.0
 153      &      .and. mod(itime,count_reset_vel).eq.0) then
 154           call random_vel
 155           write(iout,'(a,f20.2)')
 156      &     "Velocities reset to random values, time",totT
 157           do i=0,2*nres
 158             do j=1,3
 159               d_t_old(j,i)=d_t(j,i)
 160             enddo
 161           enddo
 162         endif
 163         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 164           call inertia_tensor
 165           call vcm_vel(vcm)
 166           do j=1,3
 167              d_t(j,0)=d_t(j,0)-vcm(j)
 168           enddo
 169           call kinetic(EK)
 170           kinetic_T=2.0d0/(dimen3*Rb)*EK
 171           scalfac=dsqrt(T_bath/kinetic_T)
 172           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 173           do i=0,2*nres
 174             do j=1,3
 175               d_t_old(j,i)=scalfac*d_t(j,i)
 176             enddo
 177           enddo
 178         endif
 179         if (lang.ne.4) then
 180           if (RESPA) then
 181 c Time-reversible RESPA algorithm
 182 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 183             call RESPA_step(itime)
 184           else
 185 c Variable time step algorithm.
 186             call velverlet_step(itime)
 187           endif
 188         else
 189 #ifdef BROWN
 190           call brown_step(itime)
 191 #else
 192           print *,"Brown dynamics not here!"
 193 #ifdef MPI
 194           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 195 #endif
 196           stop
 197 #endif
 198         endif
 199         if (ntwe.ne.0) then
 200          if (mod(itime,ntwe).eq.0) then
 201            call statout(itime)
 202 C           call enerprint(potEcomp)
 203 C           print *,itime,'AFM',Eafmforc,etot
 204          endif
 205 #ifdef VOUT
 206         do j=1,3
 207           v_work(j)=d_t(j,0)
 208         enddo
 209         ind=3
 210         do i=nnt,nct-1
 211           do j=1,3
 212             ind=ind+1
 213             v_work(ind)=d_t(j,i)
 214           enddo
 215         enddo
 216         do i=nnt,nct
 217           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 218             do j=1,3
 219               ind=ind+1
 220               v_work(ind)=d_t(j,i+nres)
 221             enddo
 222           endif
 223         enddo
 224
 225         write (66,'(80f10.5)')
 226      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 227         do i=1,ind
 228           v_transf(i)=0.0d0
 229           do j=1,ind
 230             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 231           enddo
 232            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 233         enddo
 234         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 235 #endif
 236         endif
 237         if (mod(itime,ntwx).eq.0) then
 238 c          write(iout,*) 'time=',itime
 239 C          call check_ecartint
 240           call returnbox
 241           write (tytul,'("time",f8.2)') totT
 242           if(mdpdb) then
 243              call hairpin(.true.,nharp,iharp)
 244              call secondary2(.true.)
 245              call pdbout(potE,tytul,ipdb)
 246           else
 247              call cartout(totT)
 248           endif
 249         endif
 250         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
 251            open(irest2,file=rest2name,status='unknown')
 252            write(irest2,*) totT,EK,potE,totE,t_bath
 253            do i=1,2*nres
 254             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 255            enddo
 256            do i=1,2*nres
 257             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 258            enddo
 259           close(irest2)
 260           rstcount=0
 261         endif
 262       enddo
 263 #ifdef MPI
 264       t_MD=MPI_Wtime()-tt0
 265 #else
 266       t_MD=tcpu()-tt0
 267 #endif
 268       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 269      &  '  Timing  ',
 270      & 'MD calculations setup:',t_MDsetup,
 271      & 'Energy & gradient evaluation:',t_enegrad,
 272      & 'Stochastic MD setup:',t_langsetup,
 273      & 'Stochastic MD step setup:',t_sdsetup,
 274      & 'MD steps:',t_MD
 275       write (iout,'(/28(1h=),a25,27(1h=))')
 276      & '  End of MD calculation  '
 277 #ifdef TIMING_ENE
 278       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
 279      &  " eshort",t_eshort
 280       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
 281      & " time_fricmatmult",time_fricmatmult," time_fsample ",
 282      & time_fsample
 283 #endif
 284       return
 285       end
 286 c-------------------------------------------------------------------------------
 287       subroutine velverlet_step(itime)
 288 c-------------------------------------------------------------------------------
 289 c  Perform a single velocity Verlet step; the time step can be rescaled if
 290 c  increments in accelerations exceed the threshold
 291 c-------------------------------------------------------------------------------
 292       implicit real*8 (a-h,o-z)
 293       include 'DIMENSIONS'
 294 #ifdef MPI
 295       include 'mpif.h'
 296       integer ierror,ierrcode
 297 #endif
 298       include 'COMMON.SETUP'
 299       include 'COMMON.CONTROL'
 300       include 'COMMON.VAR'
 301       include 'COMMON.MD'
 302       include 'COMMON.LAGRANGE'
 303 #ifndef LANG0
 304       include 'COMMON.LANGEVIN'
 305 #else
 306       include 'COMMON.LANGEVIN.lang0'
 307 #endif
 308       include 'COMMON.CHAIN'
 309       include 'COMMON.DERIV'
 310       include 'COMMON.GEO'
 311       include 'COMMON.LOCAL'
 312       include 'COMMON.INTERACT'
 313       include 'COMMON.IOUNITS'
 314       include 'COMMON.NAMES'
 315       include 'COMMON.TIME1'
 316       include 'COMMON.MUCA'
 317       double precision vcm(3),incr(3)
 318       double precision cm(3),L(3)
 319       integer ilen,count,rstcount
 320       external ilen
 321       character*50 tytul
 322       integer maxcount_scale /20/
 323       common /gucio/ cm
 324       double precision stochforcvec(MAXRES6)
 325       common /stochcalc/ stochforcvec
 326       integer itime
 327       logical scale
 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         else
 372           call verlet1
 373         endif
 374 c Build the chain from the newly calculated coordinates
 375         call chainbuild_cart
 376         if (rattle) call rattle1
 377         if (ntwe.ne.0) then
 378         if (large.and. mod(itime,ntwe).eq.0) then
 379           write (iout,*) "Cartesian and internal coordinates: step 1"
 380           call cartprint
 381           call intout
 382           write (iout,*) "dC"
 383           do i=0,nres
 384             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
 385      &      (dc(j,i+nres),j=1,3)
 386           enddo
 387           write (iout,*) "Accelerations"
 388           do i=0,nres
 389             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 390      &      (d_a(j,i+nres),j=1,3)
 391           enddo
 392           write (iout,*) "Velocities, step 1"
 393           do i=0,nres
 394             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 395      &      (d_t(j,i+nres),j=1,3)
 396           enddo
 397         endif
 398         endif
 399 #ifdef MPI
 400         tt0 = MPI_Wtime()
 401 #else
 402         tt0 = tcpu()
 403 #endif
 404 c Calculate energy and forces
 405         call zerograd
 406         call etotal(potEcomp)
 407 ! AL 4/17/17: Reduce the steps if NaNs occurred.
 408         if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
 409           d_time=d_time/2
 410           cycle
 411         endif
 412 ! end change
 413         if (large.and. mod(itime,ntwe).eq.0)
 414      &    call enerprint(potEcomp)
 415 #ifdef TIMING_ENE
 416 #ifdef MPI
 417         t_etotal=t_etotal+MPI_Wtime()-tt0
 418 #else
 419         t_etotal=t_etotal+tcpu()-tt0
 420 #endif
 421 #endif
 422         potE=potEcomp(0)-potEcomp(20)
 423         call cartgrad
 424 c Get the new accelerations
 425         call lagrangian
 426 #ifdef MPI
 427         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 428 #else
 429         t_enegrad=t_enegrad+tcpu()-tt0
 430 #endif
 431 c Determine maximum acceleration and scale down the timestep if needed
 432         call max_accel
 433         amax=amax/(itime_scal+1)**2
 434         call predict_edrift(epdrift)
 435         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
 436 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 437           scale=.true.
 438           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 439      &      /dlog(2.0d0)+1
 440           itime_scal=itime_scal+ifac_time
 441 c          fac_time=dmin1(damax/amax,0.5d0)
 442           fac_time=0.5d0**ifac_time
 443           d_time=d_time*fac_time
 444           if (lang.eq.2 .or. lang.eq.3) then
 445 #ifndef LANG0
 446 c            write (iout,*) "Calling sd_verlet_setup: 1"
 447 c Rescale the stochastic forces and recalculate or restore
 448 c the matrices of tinker integrator
 449             if (itime_scal.gt.maxflag_stoch) then
 450               if (large) write (iout,'(a,i5,a)')
 451      &         "Calculate matrices for stochastic step;",
 452      &         " itime_scal ",itime_scal
 453               if (lang.eq.2) then
 454                 call sd_verlet_p_setup
 455               else
 456                 call sd_verlet_ciccotti_setup
 457               endif
 458               write (iout,'(2a,i3,a,i3,1h.)')
 459      &         "Warning: cannot store matrices for stochastic",
 460      &         " integration because the index",itime_scal,
 461      &         " is greater than",maxflag_stoch
 462               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 463      &         " integration Langevin algorithm for better efficiency."
 464             else if (flag_stoch(itime_scal)) then
 465               if (large) write (iout,'(a,i5,a,l1)')
 466      &         "Restore matrices for stochastic step; itime_scal ",
 467      &         itime_scal," flag ",flag_stoch(itime_scal)
 468               do i=1,dimen
 469                 do j=1,dimen
 470                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 471                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 472                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 473                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 474                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 475                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 476                 enddo
 477               enddo
 478             else
 479               if (large) write (iout,'(2a,i5,a,l1)')
 480      &         "Calculate & store matrices for stochastic step;",
 481      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 482               if (lang.eq.2) then
 483                 call sd_verlet_p_setup
 484               else
 485                 call sd_verlet_ciccotti_setup
 486               endif
 487               flag_stoch(ifac_time)=.true.
 488               do i=1,dimen
 489                 do j=1,dimen
 490                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 491                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 492                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 493                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 494                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 495                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 496                 enddo
 497               enddo
 498             endif
 499             fac_time=1.0d0/dsqrt(fac_time)
 500             do i=1,dimen
 501               stochforcvec(i)=fac_time*stochforcvec(i)
 502             enddo
 503 #endif
 504           else if (lang.eq.1) then
 505 c Rescale the accelerations due to stochastic forces
 506             fac_time=1.0d0/dsqrt(fac_time)
 507             do i=1,dimen
 508               d_as_work(i)=d_as_work(i)*fac_time
 509             enddo
 510           endif
 511           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 512      &      "itime",itime," Timestep scaled down to ",
 513      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 514         else
 515 c Second step of the velocity Verlet algorithm
 516           if (lang.eq.2) then
 517 #ifndef LANG0
 518             call sd_verlet2
 519 #endif
 520           else if (lang.eq.3) then
 521 #ifndef LANG0
 522             call sd_verlet2_ciccotti
 523 #endif
 524           else if (lang.eq.1) then
 525             call sddir_verlet2
 526           else
 527             call verlet2
 528           endif
 529           if (rattle) call rattle2
 530           totT=totT+d_time
 531           totTafm=totT
 532 C          print *,totTafm,"TU?"
 533           if (d_time.ne.d_time0) then
 534             d_time=d_time0
 535 #ifndef   LANG0
 536             if (lang.eq.2 .or. lang.eq.3) then
 537               if (large) write (iout,'(a)')
 538      &         "Restore original matrices for stochastic step"
 539 c              write (iout,*) "Calling sd_verlet_setup: 2"
 540 c Restore the matrices of tinker integrator if the time step has been restored
 541               do i=1,dimen
 542                 do j=1,dimen
 543                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 544                   afric_mat(i,j)=afric0_mat(i,j,0)
 545                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 546                   prand_mat(i,j)=prand0_mat(i,j,0)
 547                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 548                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 549                 enddo
 550               enddo
 551             endif
 552 #endif
 553           endif
 554           scale=.false.
 555         endif
 556       enddo
 557 c Calculate the kinetic and the total energy and the kinetic temperature
 558       call kinetic(EK)
 559       totE=EK+potE
 560 c diagnostics
 561 c      call kinetic1(EK1)
 562 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 563 c end diagnostics
 564 c Couple the system to Berendsen bath if needed
 565       if (tbf .and. lang.eq.0) then
 566         call verlet_bath
 567       endif
 568       kinetic_T=2.0d0/(dimen3*Rb)*EK
 569 c Backup the coordinates, velocities, and accelerations
 570       do i=0,2*nres
 571         do j=1,3
 572           dc_old(j,i)=dc(j,i)
 573           d_t_old(j,i)=d_t(j,i)
 574           d_a_old(j,i)=d_a(j,i)
 575         enddo
 576       enddo
 577       if (ntwe.ne.0) then
 578       if (mod(itime,ntwe).eq.0 .and. large) then
 579         write (iout,*) "Velocities, step 2"
 580         do i=0,nres
 581           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 582      &      (d_t(j,i+nres),j=1,3)
 583         enddo
 584       endif
 585       endif
 586       return
 587       end
 588 c-------------------------------------------------------------------------------
 589       subroutine RESPA_step(itime)
 590 c-------------------------------------------------------------------------------
 591 c  Perform a single RESPA step.
 592 c-------------------------------------------------------------------------------
 593       implicit real*8 (a-h,o-z)
 594       include 'DIMENSIONS'
 595 #ifdef MPI
 596       include 'mpif.h'
 597       integer IERROR,ERRCODE
 598 #endif
 599       include 'COMMON.SETUP'
 600       include 'COMMON.CONTROL'
 601       include 'COMMON.VAR'
 602       include 'COMMON.MD'
 603       include 'COMMON.LAGRANGE'
 604 #ifndef LANG0
 605       include 'COMMON.LANGEVIN'
 606 #else
 607       include 'COMMON.LANGEVIN.lang0'
 608 #endif
 609       include 'COMMON.CHAIN'
 610       include 'COMMON.DERIV'
 611       include 'COMMON.GEO'
 612       include 'COMMON.LOCAL'
 613       include 'COMMON.INTERACT'
 614       include 'COMMON.IOUNITS'
 615       include 'COMMON.NAMES'
 616       include 'COMMON.TIME1'
 617       logical lprint_short
 618       common /shortcheck/ lprint_short
 619       double precision energia_short(0:n_ene),
 620      & energia_long(0:n_ene)
 621       double precision cm(3),L(3),vcm(3),incr(3)
 622       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
 623      & d_a_old0(3,0:maxres2)
 624       logical PRINT_AMTS_MSG /.false./
 625       integer ilen,count,rstcount
 626       external ilen
 627       character*50 tytul
 628       integer maxcount_scale /10/
 629       common /gucio/ cm
 630       double precision stochforcvec(MAXRES6)
 631       common /stochcalc/ stochforcvec
 632       integer itime
 633       logical scale
 634       common /cipiszcze/ itt
 635       itt=itime
 636       if (ntwe.ne.0) then
 637       if (large.and. mod(itime,ntwe).eq.0) then
 638         write (iout,*) "***************** RESPA itime",itime
 639         write (iout,*) "Cartesian and internal coordinates: step 0"
 640 c        call cartprint
 641         call pdbout(0.0d0,
 642      &  "cipiszcze                                         ",iout)
 643         call intout
 644         write (iout,*) "Accelerations from long-range forces"
 645         do i=0,nres
 646           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 647      &      (d_a(j,i+nres),j=1,3)
 648         enddo
 649         write (iout,*) "Velocities, step 0"
 650         do i=0,nres
 651           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 652      &      (d_t(j,i+nres),j=1,3)
 653         enddo
 654       endif
 655       endif
 656 c
 657 c Perform the initial RESPA step (increment velocities)
 658 c      write (iout,*) "*********************** RESPA ini"
 659       call RESPA_vel
 660       if (ntwe.ne.0) then
 661       if (mod(itime,ntwe).eq.0 .and. large) then
 662         write (iout,*) "Velocities, end"
 663         do i=0,nres
 664           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 665      &      (d_t(j,i+nres),j=1,3)
 666         enddo
 667       endif
 668       endif
 669 c Compute the short-range forces
 670 #ifdef MPI
 671       tt0 =MPI_Wtime()
 672 #else
 673       tt0 = tcpu()
 674 #endif
 675 C 7/2/2009 commented out
 676 c      call zerograd
 677 c      call etotal_short(energia_short)
 678 c      call cartgrad
 679 c      call lagrangian
 680 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
 681         do i=0,2*nres
 682           do j=1,3
 683             d_a(j,i)=d_a_short(j,i)
 684           enddo
 685         enddo
 686       if (ntwe.ne.0) then
 687       if (large.and. mod(itime,ntwe).eq.0) then
 688         write (iout,*) "energia_short",energia_short(0)
 689         write (iout,*) "Accelerations from short-range forces"
 690         do i=0,nres
 691           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 692      &      (d_a(j,i+nres),j=1,3)
 693         enddo
 694       endif
 695       endif
 696 #ifdef MPI
 697         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 698 #else
 699         t_enegrad=t_enegrad+tcpu()-tt0
 700 #endif
 701       do i=0,2*nres
 702         do j=1,3
 703           dc_old(j,i)=dc(j,i)
 704           d_t_old(j,i)=d_t(j,i)
 705           d_a_old(j,i)=d_a(j,i)
 706         enddo
 707       enddo
 708 c 6/30/08 A-MTS: attempt at increasing the split number
 709       do i=0,2*nres
 710         do j=1,3
 711           dc_old0(j,i)=dc_old(j,i)
 712           d_t_old0(j,i)=d_t_old(j,i)
 713           d_a_old0(j,i)=d_a_old(j,i)
 714         enddo
 715       enddo
 716       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 717       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 718 c
 719       scale=.true.
 720       d_time0=d_time
 721       do while (scale)
 722
 723       scale=.false.
 724 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 725 c Split the time step
 726       d_time=d_time0/ntime_split
 727 c Perform the short-range RESPA steps (velocity Verlet increments of
 728 c positions and velocities using short-range forces)
 729 c      write (iout,*) "*********************** RESPA split"
 730       do itsplit=1,ntime_split
 731         if (lang.eq.1) then
 732           call sddir_precalc
 733         else if (lang.eq.2 .or. lang.eq.3) then
 734 #ifndef LANG0
 735           call stochastic_force(stochforcvec)
 736 #else
 737           write (iout,*)
 738      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 739 #ifdef MPI
 740           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 741 #endif
 742           stop
 743 #endif
 744         endif
 745 c First step of the velocity Verlet algorithm
 746         if (lang.eq.2) then
 747 #ifndef LANG0
 748           call sd_verlet1
 749 #endif
 750         else if (lang.eq.3) then
 751 #ifndef LANG0
 752           call sd_verlet1_ciccotti
 753 #endif
 754         else if (lang.eq.1) then
 755           call sddir_verlet1
 756         else
 757           call verlet1
 758         endif
 759 c Build the chain from the newly calculated coordinates
 760         call chainbuild_cart
 761         if (rattle) call rattle1
 762         if (ntwe.ne.0) then
 763         if (large.and. mod(itime,ntwe).eq.0) then
 764           write (iout,*) "***** ITSPLIT",itsplit
 765           write (iout,*) "Cartesian and internal coordinates: step 1"
 766           call pdbout(0.0d0,
 767      &     "cipiszcze                                         ",iout)
 768 c          call cartprint
 769           call intout
 770           write (iout,*) "Velocities, step 1"
 771           do i=0,nres
 772             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 773      &        (d_t(j,i+nres),j=1,3)
 774           enddo
 775         endif
 776         endif
 777 #ifdef MPI
 778         tt0 = MPI_Wtime()
 779 #else
 780         tt0 = tcpu()
 781 #endif
 782 c Calculate energy and forces
 783 c        if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
 784         call zerograd
 785         call etotal_short(energia_short)
 786 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
 787         if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
 788      &  then
 789           if (PRINT_AMTS_MSG) write (iout,*)
 790      &     "Infinities/NaNs in energia_short",
 791      &      energia_short(0),"; increasing ntime_split to",ntime_split
 792           ntime_split=ntime_split*2
 793           if (ntime_split.gt.maxtime_split) then
 794 #ifdef MPI
 795           write (iout,*) "Cannot rescue the run; aborting job.",
 796      &      " Retry with a smaller time step"
 797           call flush(iout)
 798           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 799 #else
 800           write (iout,*) "Cannot rescue the run; terminating.",
 801      &      " Retry with a smaller time step"
 802 #endif
 803           endif
 804           exit
 805         endif
 806 ! End change
 807         if (large.and. mod(itime,ntwe).eq.0)
 808      &    call enerprint(energia_short)
 809 #ifdef TIMING_ENE
 810 #ifdef MPI
 811         t_eshort=t_eshort+MPI_Wtime()-tt0
 812 #else
 813         t_eshort=t_eshort+tcpu()-tt0
 814 #endif
 815 #endif
 816         call cartgrad
 817 c Get the new accelerations
 818         call lagrangian
 819 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 820         do i=0,2*nres
 821           do j=1,3
 822             d_a_short(j,i)=d_a(j,i)
 823           enddo
 824         enddo
 825         if (ntwe.ne.0) then
 826         if (large.and. mod(itime,ntwe).eq.0) then
 827           write (iout,*)"energia_short",energia_short(0)
 828           write (iout,*) "Accelerations from short-range forces"
 829           do i=0,nres
 830             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 831      &        (d_a(j,i+nres),j=1,3)
 832           enddo
 833         endif
 834         endif
 835 c 6/30/08 A-MTS
 836 c Determine maximum acceleration and scale down the timestep if needed
 837         call max_accel
 838         amax=amax/ntime_split**2
 839         call predict_edrift(epdrift)
 840         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 841      &   write (iout,*) "amax",amax," damax",damax,
 842      &   " epdrift",epdrift," epdriftmax",epdriftmax
 843 c Exit loop and try with increased split number if the change of
 844 c acceleration is too big
 845         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 846           if (ntime_split.lt.maxtime_split) then
 847             scale=.true.
 848             ntime_split=ntime_split*2
 849             do i=0,2*nres
 850               do j=1,3
 851                 dc_old(j,i)=dc_old0(j,i)
 852                 d_t_old(j,i)=d_t_old0(j,i)
 853                 d_a_old(j,i)=d_a_old0(j,i)
 854               enddo
 855             enddo
 856             if (PRINT_AMTS_MSG) then
 857             write (iout,*) "acceleration/energy drift too large",amax,
 858      &      epdrift," split increased to ",ntime_split," itime",itime,
 859      &       " itsplit",itsplit
 860             endif
 861             exit
 862           else
 863             write (iout,*)
 864      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 865      &       maxtime_split,
 866      &      " already reached!!! Trying to carry on!"
 867           endif
 868         endif
 869 #ifdef MPI
 870         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 871 #else
 872         t_enegrad=t_enegrad+tcpu()-tt0
 873 #endif
 874 c Second step of the velocity Verlet algorithm
 875         if (lang.eq.2) then
 876 #ifndef LANG0
 877           call sd_verlet2
 878 #endif
 879         else if (lang.eq.3) then
 880 #ifndef LANG0
 881           call sd_verlet2_ciccotti
 882 #endif
 883         else if (lang.eq.1) then
 884           call sddir_verlet2
 885         else
 886           call verlet2
 887         endif
 888         if (rattle) call rattle2
 889 c Backup the coordinates, velocities, and accelerations
 890         do i=0,2*nres
 891           do j=1,3
 892             dc_old(j,i)=dc(j,i)
 893             d_t_old(j,i)=d_t(j,i)
 894             d_a_old(j,i)=d_a(j,i)
 895           enddo
 896         enddo
 897       enddo
 898
 899       enddo ! while scale
 900
 901 c Restore the time step
 902       d_time=d_time0
 903 c Compute long-range forces
 904 #ifdef MPI
 905       tt0 =MPI_Wtime()
 906 #else
 907       tt0 = tcpu()
 908 #endif
 909       call zerograd
 910       call etotal_long(energia_long)
 911 ! AL 4/17/2017 Handling NaNs
 912       if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
 913 #ifdef MPI
 914         write (iout,*)
 915      &   "Infinitied/NaNs in energia_long, Aborting MPI job."
 916         call enerprint(energia_long)
 917         call flush(iout)
 918         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 919 #else
 920         write (iout,*) "Infinitied/NaNs in energia_long, terminating."
 921         call enerprint(energia_long)
 922         stop
 923 #endif
 924       endif
 925 ! end change
 926 c        lprint_short=.false.
 927       if (large.and. mod(itime,ntwe).eq.0)
 928      &    call enerprint(energia_long)
 929 #ifdef TIMING_ENE
 930 #ifdef MPI
 931         t_elong=t_elong+MPI_Wtime()-tt0
 932 #else
 933         t_elong=t_elong+tcpu()-tt0
 934 #endif
 935 #endif
 936       call cartgrad
 937       call lagrangian
 938 #ifdef MPI
 939         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 940 #else
 941         t_enegrad=t_enegrad+tcpu()-tt0
 942 #endif
 943 c Compute accelerations from long-range forces
 944       if (ntwe.ne.0) then
 945       if (large.and. mod(itime,ntwe).eq.0) then
 946         write (iout,*) "energia_long",energia_long(0)
 947         write (iout,*) "Cartesian and internal coordinates: step 2"
 948 c        call cartprint
 949         call pdbout(0.0d0,
 950      &    cipiszcze                                         ,iout)
 951         call intout
 952         write (iout,*) "Accelerations from long-range forces"
 953         do i=0,nres
 954           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 955      &      (d_a(j,i+nres),j=1,3)
 956         enddo
 957         write (iout,*) "Velocities, step 2"
 958         do i=0,nres
 959           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 960      &      (d_t(j,i+nres),j=1,3)
 961         enddo
 962       endif
 963       endif
 964 c Compute the final RESPA step (increment velocities)
 965 c      write (iout,*) "*********************** RESPA fin"
 966       call RESPA_vel
 967 c Compute the complete potential energy
 968       do i=0,n_ene
 969         potEcomp(i)=energia_short(i)+energia_long(i)
 970       enddo
 971       potE=potEcomp(0)-potEcomp(20)
 972       if (ntwe.ne.0) then
 973       if (large.and. mod(itime,ntwe).eq.0) then
 974         call enerprint(potEcomp)
 975         write (iout,*) "potE",potD
 976       endif
 977       endif
 978 c      potE=energia_short(0)+energia_long(0)
 979       totT=totT+d_time
 980       totTafm=totT
 981 c Calculate the kinetic and the total energy and the kinetic temperature
 982       call kinetic(EK)
 983       totE=EK+potE
 984 c Couple the system to Berendsen bath if needed
 985       if (tbf .and. lang.eq.0) then
 986         call verlet_bath
 987       endif
 988       kinetic_T=2.0d0/(dimen3*Rb)*EK
 989 c Backup the coordinates, velocities, and accelerations
 990       if (ntwe.ne.0) then
 991       if (mod(itime,ntwe).eq.0 .and. large) then
 992         write (iout,*) "Velocities, end"
 993         do i=0,nres
 994           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 995      &      (d_t(j,i+nres),j=1,3)
 996         enddo
 997       endif
 998       endif
 999       return
1000       end
1001 c---------------------------------------------------------------------
1002       subroutine RESPA_vel
1003 c  First and last RESPA step (incrementing velocities using long-range
1004 c  forces).
1005       implicit real*8 (a-h,o-z)
1006       include 'DIMENSIONS'
1007       include 'COMMON.CONTROL'
1008       include 'COMMON.VAR'
1009       include 'COMMON.MD'
1010       include 'COMMON.LAGRANGE'
1011       include 'COMMON.CHAIN'
1012       include 'COMMON.DERIV'
1013       include 'COMMON.GEO'
1014       include 'COMMON.LOCAL'
1015       include 'COMMON.INTERACT'
1016       include 'COMMON.IOUNITS'
1017       include 'COMMON.NAMES'
1018       do j=1,3
1019         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1020       enddo
1021       do i=nnt,nct-1
1022         do j=1,3
1023           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1024         enddo
1025       enddo
1026       do i=nnt,nct
1027         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1028           inres=i+nres
1029           do j=1,3
1030             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1031           enddo
1032         endif
1033       enddo
1034       return
1035       end
1036 c-----------------------------------------------------------------
1037       subroutine verlet1
1038 c Applying velocity Verlet algorithm - step 1 to coordinates
1039       implicit real*8 (a-h,o-z)
1040       include 'DIMENSIONS'
1041       include 'COMMON.CONTROL'
1042       include 'COMMON.VAR'
1043       include 'COMMON.MD'
1044       include 'COMMON.LAGRANGE'
1045       include 'COMMON.CHAIN'
1046       include 'COMMON.DERIV'
1047       include 'COMMON.GEO'
1048       include 'COMMON.LOCAL'
1049       include 'COMMON.INTERACT'
1050       include 'COMMON.IOUNITS'
1051       include 'COMMON.NAMES'
1052       double precision adt,adt2
1053
1054 #ifdef DEBUG
1055       write (iout,*) "VELVERLET1 START: DC"
1056       do i=0,nres
1057         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1058      &   (dc(j,i+nres),j=1,3)
1059       enddo
1060 #endif
1061       do j=1,3
1062         adt=d_a_old(j,0)*d_time
1063         adt2=0.5d0*adt
1064         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1065         d_t_new(j,0)=d_t_old(j,0)+adt2
1066         d_t(j,0)=d_t_old(j,0)+adt
1067       enddo
1068       do i=nnt,nct-1
1069 C      SPYTAC ADAMA
1070 C       do i=0,nres
1071         do j=1,3
1072           adt=d_a_old(j,i)*d_time
1073           adt2=0.5d0*adt
1074           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1075           d_t_new(j,i)=d_t_old(j,i)+adt2
1076           d_t(j,i)=d_t_old(j,i)+adt
1077         enddo
1078       enddo
1079       do i=nnt,nct
1080 C        do i=0,nres
1081         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1082           inres=i+nres
1083           do j=1,3
1084             adt=d_a_old(j,inres)*d_time
1085             adt2=0.5d0*adt
1086             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1087             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1088             d_t(j,inres)=d_t_old(j,inres)+adt
1089           enddo
1090         endif
1091       enddo
1092 #ifdef DEBUG
1093       write (iout,*) "VELVERLET1 END: DC"
1094       do i=0,nres
1095         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1096      &   (dc(j,i+nres),j=1,3)
1097       enddo
1098 #endif
1099       return
1100       end
1101 c---------------------------------------------------------------------
1102       subroutine verlet2
1103 c  Step 2 of the velocity Verlet algorithm: update velocities
1104       implicit real*8 (a-h,o-z)
1105       include 'DIMENSIONS'
1106       include 'COMMON.CONTROL'
1107       include 'COMMON.VAR'
1108       include 'COMMON.MD'
1109       include 'COMMON.LAGRANGE'
1110       include 'COMMON.CHAIN'
1111       include 'COMMON.DERIV'
1112       include 'COMMON.GEO'
1113       include 'COMMON.LOCAL'
1114       include 'COMMON.INTERACT'
1115       include 'COMMON.IOUNITS'
1116       include 'COMMON.NAMES'
1117       do j=1,3
1118         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1119       enddo
1120       do i=nnt,nct-1
1121         do j=1,3
1122           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1123         enddo
1124       enddo
1125       do i=nnt,nct
1126         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1127           inres=i+nres
1128           do j=1,3
1129             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1130           enddo
1131         endif
1132       enddo
1133       return
1134       end
1135 c-----------------------------------------------------------------
1136       subroutine sddir_precalc
1137 c Applying velocity Verlet algorithm - step 1 to coordinates
1138       implicit real*8 (a-h,o-z)
1139       include 'DIMENSIONS'
1140 #ifdef MPI
1141       include 'mpif.h'
1142 #endif
1143       include 'COMMON.CONTROL'
1144       include 'COMMON.VAR'
1145       include 'COMMON.MD'
1146       include 'COMMON.LAGRANGE'
1147 #ifndef LANG0
1148       include 'COMMON.LANGEVIN'
1149 #else
1150       include 'COMMON.LANGEVIN.lang0'
1151 #endif
1152       include 'COMMON.CHAIN'
1153       include 'COMMON.DERIV'
1154       include 'COMMON.GEO'
1155       include 'COMMON.LOCAL'
1156       include 'COMMON.INTERACT'
1157       include 'COMMON.IOUNITS'
1158       include 'COMMON.NAMES'
1159       include 'COMMON.TIME1'
1160       double precision stochforcvec(MAXRES6)
1161       common /stochcalc/ stochforcvec
1162 c
1163 c Compute friction and stochastic forces
1164 c
1165 #ifdef MPI
1166       time00=MPI_Wtime()
1167 #else
1168       time00=tcpu()
1169 #endif
1170       call friction_force
1171 #ifdef MPI
1172       time_fric=time_fric+MPI_Wtime()-time00
1173       time00=MPI_Wtime()
1174 #else
1175       time_fric=time_fric+tcpu()-time00
1176       time00=tcpu()
1177 #endif
1178       call stochastic_force(stochforcvec)
1179 #ifdef MPI
1180       time_stoch=time_stoch+MPI_Wtime()-time00
1181 #else
1182       time_stoch=time_stoch+tcpu()-time00
1183 #endif
1184 c
1185 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1186 c forces (d_as_work)
1187 c
1188       call ginv_mult(fric_work, d_af_work)
1189       call ginv_mult(stochforcvec, d_as_work)
1190       return
1191       end
1192 c---------------------------------------------------------------------
1193       subroutine sddir_verlet1
1194 c Applying velocity Verlet algorithm - step 1 to velocities
1195       implicit real*8 (a-h,o-z)
1196       include 'DIMENSIONS'
1197       include 'COMMON.CONTROL'
1198       include 'COMMON.VAR'
1199       include 'COMMON.MD'
1200       include 'COMMON.LAGRANGE'
1201 #ifndef LANG0
1202       include 'COMMON.LANGEVIN'
1203 #else
1204       include 'COMMON.LANGEVIN.lang0'
1205 #endif
1206       include 'COMMON.CHAIN'
1207       include 'COMMON.DERIV'
1208       include 'COMMON.GEO'
1209       include 'COMMON.LOCAL'
1210       include 'COMMON.INTERACT'
1211       include 'COMMON.IOUNITS'
1212       include 'COMMON.NAMES'
1213 c Revised 3/31/05 AL: correlation between random contributions to
1214 c position and velocity increments included.
1215       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1216       double precision adt,adt2
1217 c
1218 c Add the contribution from BOTH friction and stochastic force to the
1219 c coordinates, but ONLY the contribution from the friction forces to velocities
1220 c
1221       do j=1,3
1222         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1223         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1224         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1225         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1226         d_t(j,0)=d_t_old(j,0)+adt
1227       enddo
1228       ind=3
1229       do i=nnt,nct-1
1230         do j=1,3
1231           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1232           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1233           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1234           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1235           d_t(j,i)=d_t_old(j,i)+adt
1236         enddo
1237         ind=ind+3
1238       enddo
1239       do i=nnt,nct
1240         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1241           inres=i+nres
1242           do j=1,3
1243             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1244             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1245             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1246             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1247             d_t(j,inres)=d_t_old(j,inres)+adt
1248           enddo
1249           ind=ind+3
1250         endif
1251       enddo
1252       return
1253       end
1254 c---------------------------------------------------------------------
1255       subroutine sddir_verlet2
1256 c  Calculating the adjusted velocities for accelerations
1257       implicit real*8 (a-h,o-z)
1258       include 'DIMENSIONS'
1259       include 'COMMON.CONTROL'
1260       include 'COMMON.VAR'
1261       include 'COMMON.MD'
1262       include 'COMMON.LAGRANGE'
1263 #ifndef LANG0
1264       include 'COMMON.LANGEVIN'
1265 #else
1266       include 'COMMON.LANGEVIN.lang0'
1267 #endif
1268       include 'COMMON.CHAIN'
1269       include 'COMMON.DERIV'
1270       include 'COMMON.GEO'
1271       include 'COMMON.LOCAL'
1272       include 'COMMON.INTERACT'
1273       include 'COMMON.IOUNITS'
1274       include 'COMMON.NAMES'
1275       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1276       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1277 c Revised 3/31/05 AL: correlation between random contributions to
1278 c position and velocity increments included.
1279 c The correlation coefficients are calculated at low-friction limit.
1280 c Also, friction forces are now not calculated with new velocities.
1281
1282 c      call friction_force
1283       call stochastic_force(stochforcvec)
1284 c
1285 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1286 c forces (d_as_work)
1287 c
1288       call ginv_mult(stochforcvec, d_as_work1)
1289
1290 c
1291 c Update velocities
1292 c
1293       do j=1,3
1294         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1295      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1296       enddo
1297       ind=3
1298       do i=nnt,nct-1
1299         do j=1,3
1300           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1301      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1302         enddo
1303         ind=ind+3
1304       enddo
1305       do i=nnt,nct
1306         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1307           inres=i+nres
1308           do j=1,3
1309             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1310      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1311      &       +cos60*d_as_work1(ind+j))*d_time
1312           enddo
1313           ind=ind+3
1314         endif
1315       enddo
1316       return
1317       end
1318 c---------------------------------------------------------------------
1319       subroutine max_accel
1320 c
1321 c Find the maximum difference in the accelerations of the the sites
1322 c at the beginning and the end of the time step.
1323 c
1324       implicit real*8 (a-h,o-z)
1325       include 'DIMENSIONS'
1326       include 'COMMON.CONTROL'
1327       include 'COMMON.VAR'
1328       include 'COMMON.MD'
1329       include 'COMMON.LAGRANGE'
1330       include 'COMMON.CHAIN'
1331       include 'COMMON.DERIV'
1332       include 'COMMON.GEO'
1333       include 'COMMON.LOCAL'
1334       include 'COMMON.INTERACT'
1335       include 'COMMON.IOUNITS'
1336       double precision aux(3),accel(3),accel_old(3),dacc
1337       do j=1,3
1338 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1339          accel_old(j)=d_a_old(j,0)
1340          accel(j)=d_a(j,0)
1341       enddo
1342       amax=0.0d0
1343       do i=nnt,nct
1344 c Backbone
1345         if (i.lt.nct) then
1346 c 7/3/08 changed to asymmetric difference
1347           do j=1,3
1348 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1349             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1350             accel(j)=accel(j)+0.5d0*d_a(j,i)
1351 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1352             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1353               dacc=dabs(accel(j)-accel_old(j))
1354 c              write (iout,*) i,dacc
1355               if (dacc.gt.amax) amax=dacc
1356             endif
1357           enddo
1358         endif
1359       enddo
1360 c Side chains
1361       do j=1,3
1362 c        accel(j)=aux(j)
1363         accel_old(j)=d_a_old(j,0)
1364         accel(j)=d_a(j,0)
1365       enddo
1366       if (nnt.eq.2) then
1367         do j=1,3
1368           accel_old(j)=accel_old(j)+d_a_old(j,1)
1369           accel(j)=accel(j)+d_a(j,1)
1370         enddo
1371       endif
1372       do i=nnt,nct
1373         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1374           do j=1,3
1375 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1376             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1377             accel(j)=accel(j)+d_a(j,i+nres)
1378           enddo
1379         endif
1380         do j=1,3
1381 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1382           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1383             dacc=dabs(accel(j)-accel_old(j))
1384 c            write (iout,*) "side-chain",i,dacc
1385             if (dacc.gt.amax) amax=dacc
1386           endif
1387         enddo
1388         do j=1,3
1389           accel_old(j)=accel_old(j)+d_a_old(j,i)
1390           accel(j)=accel(j)+d_a(j,i)
1391 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1392         enddo
1393       enddo
1394       return
1395       end
1396 c---------------------------------------------------------------------
1397       subroutine predict_edrift(epdrift)
1398 c
1399 c Predict the drift of the potential energy
1400 c
1401       implicit real*8 (a-h,o-z)
1402       include 'DIMENSIONS'
1403       include 'COMMON.CONTROL'
1404       include 'COMMON.VAR'
1405       include 'COMMON.MD'
1406       include 'COMMON.LAGRANGE'
1407       include 'COMMON.CHAIN'
1408       include 'COMMON.DERIV'
1409       include 'COMMON.GEO'
1410       include 'COMMON.LOCAL'
1411       include 'COMMON.INTERACT'
1412       include 'COMMON.IOUNITS'
1413       include 'COMMON.MUCA'
1414       double precision epdrift,epdriftij
1415 c Drift of the potential energy
1416       epdrift=0.0d0
1417       do i=nnt,nct
1418 c Backbone
1419         if (i.lt.nct) then
1420           do j=1,3
1421             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1422             if (lmuca) epdriftij=epdriftij*factor
1423 c            write (iout,*) "back",i,j,epdriftij
1424             if (epdriftij.gt.epdrift) epdrift=epdriftij
1425           enddo
1426         endif
1427 c Side chains
1428         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1429           do j=1,3
1430             epdriftij=
1431      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1432             if (lmuca) epdriftij=epdriftij*factor
1433 c            write (iout,*) "side",i,j,epdriftij
1434             if (epdriftij.gt.epdrift) epdrift=epdriftij
1435           enddo
1436         endif
1437       enddo
1438       epdrift=0.5d0*epdrift*d_time*d_time
1439 c      write (iout,*) "epdrift",epdrift
1440       return
1441       end
1442 c-----------------------------------------------------------------------
1443       subroutine verlet_bath
1444 c
1445 c  Coupling to the thermostat by using the Berendsen algorithm
1446 c
1447       implicit real*8 (a-h,o-z)
1448       include 'DIMENSIONS'
1449       include 'COMMON.CONTROL'
1450       include 'COMMON.VAR'
1451       include 'COMMON.MD'
1452       include 'COMMON.LAGRANGE'
1453       include 'COMMON.CHAIN'
1454       include 'COMMON.DERIV'
1455       include 'COMMON.GEO'
1456       include 'COMMON.LOCAL'
1457       include 'COMMON.INTERACT'
1458       include 'COMMON.IOUNITS'
1459       include 'COMMON.NAMES'
1460       double precision T_half,fact
1461 c
1462       T_half=2.0d0/(dimen3*Rb)*EK
1463       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1464 c      write(iout,*) "T_half", T_half
1465 c      write(iout,*) "EK", EK
1466 c      write(iout,*) "fact", fact
1467       do j=1,3
1468         d_t(j,0)=fact*d_t(j,0)
1469       enddo
1470       do i=nnt,nct-1
1471         do j=1,3
1472           d_t(j,i)=fact*d_t(j,i)
1473         enddo
1474       enddo
1475       do i=nnt,nct
1476         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1477           inres=i+nres
1478           do j=1,3
1479             d_t(j,inres)=fact*d_t(j,inres)
1480           enddo
1481         endif
1482       enddo
1483       return
1484       end
1485 c---------------------------------------------------------
1486       subroutine init_MD
1487 c  Set up the initial conditions of a MD simulation
1488       implicit real*8 (a-h,o-z)
1489       include 'DIMENSIONS'
1490 #ifdef MP
1491       include 'mpif.h'
1492       character*16 form
1493       integer IERROR,ERRCODE
1494 #endif
1495       include 'COMMON.SETUP'
1496       include 'COMMON.CONTROL'
1497       include 'COMMON.VAR'
1498       include 'COMMON.MD'
1499       include 'COMMON.QRESTR'
1500       include 'COMMON.LAGRANGE'
1501 #ifndef LANG0
1502       include 'COMMON.LANGEVIN'
1503 #else
1504       include 'COMMON.LANGEVIN.lang0'
1505 #endif
1506       include 'COMMON.CHAIN'
1507       include 'COMMON.DERIV'
1508       include 'COMMON.GEO'
1509       include 'COMMON.LOCAL'
1510       include 'COMMON.INTERACT'
1511       include 'COMMON.IOUNITS'
1512       include 'COMMON.NAMES'
1513       include 'COMMON.REMD'
1514       real*8 energia_long(0:n_ene),
1515      &  energia_short(0:n_ene),vcm(3),incr(3)
1516       double precision cm(3),L(3),xv,sigv,lowb,highb
1517       double precision varia(maxvar),energia(0:n_ene)
1518       character*256 qstr
1519       integer ilen
1520       external ilen
1521       character*50 tytul
1522       logical file_exist
1523       common /gucio/ cm
1524       write (iout,*) "init_MD INDPDB",indpdb
1525       d_time0=d_time
1526 c      write(iout,*) "d_time", d_time
1527 c Compute the standard deviations of stochastic forces for Langevin dynamics
1528 c if the friction coefficients do not depend on surface area
1529       if (lang.gt.0 .and. .not.surfarea) then
1530         do i=nnt,nct-1
1531           stdforcp(i)=stdfp*dsqrt(gamp)
1532         enddo
1533         do i=nnt,nct
1534           if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1535      &                *dsqrt(gamsc(iabs(itype(i))))
1536         enddo
1537       endif
1538 c Open the pdb file for snapshotshots
1539 #ifdef MPI
1540       if(mdpdb) then
1541         if (ilen(tmpdir).gt.0)
1542      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1543      &      liczba(:ilen(liczba))//".pdb")
1544         open(ipdb,
1545      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1546      &  //".pdb")
1547       else
1548 #ifdef NOXDR
1549         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1550      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1551      &      liczba(:ilen(liczba))//".x")
1552         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1553      &  //".x"
1554 #else
1555         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1556      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1557      &      liczba(:ilen(liczba))//".cx")
1558         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1559      &  //".cx"
1560 #endif
1561       endif
1562 #else
1563       if(mdpdb) then
1564          if (ilen(tmpdir).gt.0)
1565      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1566          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1567       else
1568          if (ilen(tmpdir).gt.0)
1569      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1570          cartname=prefix(:ilen(prefix))//"_MD.cx"
1571       endif
1572 #endif
1573       if (usampl) then
1574         write (qstr,'(256(1h ))')
1575         ipos=1
1576         do i=1,nfrag
1577           iq = qinfrag(i,iset)*10
1578           iw = wfrag(i,iset)/100
1579           if (iw.gt.0) then
1580             if(me.eq.king.or..not.out1file)
1581      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1582             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1583             ipos=ipos+7
1584           endif
1585         enddo
1586         do i=1,npair
1587           iq = qinpair(i,iset)*10
1588           iw = wpair(i,iset)/100
1589           if (iw.gt.0) then
1590             if(me.eq.king.or..not.out1file)
1591      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1592             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1593             ipos=ipos+7
1594           endif
1595         enddo
1596 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1597 #ifdef NOXDR
1598 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1599 #else
1600 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1601 #endif
1602 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1603       endif
1604       icg=1
1605       write (iout,*) "REST ",rest
1606       if (rest) then
1607        if (restart1file) then
1608          if (me.eq.king)
1609      &     inquire(file=mremd_rst_name,exist=file_exist)
1610 #ifdef MPI
1611            write (*,*) me," Before broadcast: file_exist",file_exist
1612          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1613      &          IERR)
1614          write (*,*) me," After broadcast: file_exist",file_exist
1615 c        inquire(file=mremd_rst_name,exist=file_exist)
1616 #endif
1617         if(me.eq.king.or..not.out1file)
1618      &   write(iout,*) "Initial state read by master and distributed"
1619        else
1620          if (ilen(tmpdir).gt.0)
1621      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1622      &      //liczba(:ilen(liczba))//'.rst')
1623         inquire(file=rest2name,exist=file_exist)
1624        endif
1625        if(file_exist) then
1626          if(.not.restart1file) then
1627            if(me.eq.king.or..not.out1file)
1628      &      write(iout,*) "Initial state will be read from file ",
1629      &      rest2name(:ilen(rest2name))
1630            call readrst
1631          endif
1632          call rescale_weights(t_bath)
1633        else
1634         rest=.false.
1635         if(me.eq.king.or..not.out1file)then
1636          if (restart1file) then
1637           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1638      &       " does not exist"
1639          else
1640           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1641      &       " does not exist"
1642          endif
1643          write(iout,*) "Initial velocities randomly generated"
1644         endif
1645         call random_vel
1646         totT=0.0d0
1647         totTafm=totT
1648        endif
1649       else
1650 c Generate initial velocities
1651         if(me.eq.king.or..not.out1file)
1652      &   write(iout,*) "Initial velocities randomly generated"
1653         call random_vel
1654         totT=0.0d0
1655 CtotTafm is the variable for AFM time which eclipsed during
1656         totTafm=totT
1657       endif
1658 c      rest2name = prefix(:ilen(prefix))//'.rst'
1659       if(me.eq.king.or..not.out1file)then
1660        write (iout,*) "Initial velocities"
1661        do i=0,nres
1662          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1663      &   (d_t(j,i+nres),j=1,3)
1664        enddo
1665 c  Zeroing the total angular momentum of the system
1666        write(iout,*) "Calling the zero-angular
1667      & momentum subroutine"
1668       endif
1669       call inertia_tensor
1670 c  Getting the potential energy and forces and velocities and accelerations
1671       call vcm_vel(vcm)
1672       write (iout,*) "velocity of the center of the mass:"
1673       write (iout,*) (vcm(j),j=1,3)
1674       call flush(iout)
1675       do j=1,3
1676         d_t(j,0)=d_t(j,0)-vcm(j)
1677       enddo
1678 c Removing the velocity of the center of mass
1679       call vcm_vel(vcm)
1680       if(me.eq.king.or..not.out1file)then
1681        write (iout,*) "vcm right after adjustment:"
1682        write (iout,*) (vcm(j),j=1,3)
1683       endif
1684       call flush(iout)
1685       write (iout,*) "init_MD before initial structure REST ",rest
1686       if (.not.rest) then
1687         if (iranconf.ne.0) then
1688 c 8/22/17 AL Loop to produce a low-energy random conformation
1689           do iranmin=1,10
1690             call chainbuild
1691             if (overlapsc) then
1692               print *, 'Calling OVERLAP_SC'
1693               call overlap_sc(fail)
1694             endif
1695
1696             if (searchsc) then
1697               call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1698               print *,'SC_move',nft_sc,etot
1699               if(me.eq.king.or..not.out1file)
1700      &          write(iout,*) 'SC_move',nft_sc,etot
1701             endif
1702
1703             if (dccart) then
1704               if (me.eq.king.or..not.out1file) write(iout,*)
1705      &         'Minimizing random structure: Calling MINIM_DC'
1706                 call minim_dc(etot,iretcode,nfun)
1707             else
1708               call geom_to_var(nvar,varia)
1709               if (me.eq.king.or..not.out1file) write (iout,*)
1710      &          'Minimizing random structure: Calling MINIMIZE.'
1711               call minimize(etot,varia,iretcode,nfun)
1712               call var_to_geom(nvar,varia)
1713             endif
1714             if (me.eq.king.or..not.out1file)
1715      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1716             if (isnan(etot) .or. etot.gt.1.0d4) then
1717                write (iout,*) "Energy too large",etot,
1718      &          " trying another random conformation"
1719               do itrial=1,100
1720                 itmp=1
1721                 call gen_rand_conf(itmp,*30)
1722                 goto 40
1723    30           write (iout,*) 'Failed to generate random conformation',
1724      &            ', itrial=',itrial
1725                 write (*,*) 'Processor:',me,
1726      &            ' Failed to generate random conformation',
1727      &            ' itrial=',itrial
1728                 call intout
1729
1730 #ifdef AIX
1731                 call flush_(iout)
1732 #else
1733                 call flush(iout)
1734 #endif
1735               enddo
1736               write (iout,'(a,i3,a)') 'Processor:',me,
1737      &        ' error in generating random conformation.'
1738               write (*,'(a,i3,a)') 'Processor:',me,
1739      &        ' error in generating random conformation.'
1740               call flush(iout)
1741 #ifdef MPI
1742               call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1743 #else
1744               stop
1745 #endif
1746    40         continue
1747             else
1748               goto 44
1749             endif
1750           enddo
1751           write (iout,'(a,i3,a)') 'Processor:',me,
1752      &        ' failed to generate a low-energy random conformation.'
1753           write (*,'(a,i3,a)') 'Processor:',me,
1754      &        ' failed to generate a low-energy random conformation.'
1755             call flush(iout)
1756 #ifdef MPI
1757             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1758 #else
1759             stop
1760 #endif
1761    44     continue
1762         else if (preminim) then
1763           if (start_from_model) then
1764             i_model=iran_num(1,constr_homology)
1765             write (iout,*) 'starting from model ',i_model
1766             do i=1,2*nres
1767               do j=1,3
1768                 c(j,i)=chomo(j,i,i_model)
1769               enddo
1770             enddo
1771             call int_from_cart(.true.,.false.)
1772             call sc_loc_geom(.false.)
1773             do i=1,nres-1
1774               do j=1,3
1775                 dc(j,i)=c(j,i+1)-c(j,i)
1776                 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1777               enddo
1778             enddo
1779             do i=2,nres-1
1780               do j=1,3
1781                 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1782                 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1783               enddo
1784             enddo
1785           endif
1786 ! Remove SC overlaps if requested
1787           if (overlapsc) then
1788             write (iout,*) 'Calling OVERLAP_SC'
1789             call overlap_sc(fail)
1790           endif
1791 ! Search for better SC rotamers if requested
1792           if (searchsc) then
1793             call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1794             print *,'SC_move',nft_sc,etot
1795             if (me.eq.king.or..not.out1file)
1796      &        write(iout,*) 'SC_move',nft_sc,etot
1797           endif
1798           call etotal(energia(0))
1799 C 8/22/17 AL Minimize initial structure
1800           if (dccart) then
1801             if (me.eq.king.or..not.out1file) write(iout,*)
1802      &      'Minimizing initial PDB structure: Calling MINIM_DC'
1803             call minim_dc(etot,iretcode,nfun)
1804           else
1805             call geom_to_var(nvar,varia)
1806             if(me.eq.king.or..not.out1file) write (iout,*)
1807      &      'Minimizing initial PDB structure: Calling MINIMIZE.'
1808             call minimize(etot,varia,iretcode,nfun)
1809             call var_to_geom(nvar,varia)
1810           endif
1811           if (me.eq.king.or..not.out1file)
1812      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1813           if(me.eq.king.or..not.out1file)
1814      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1815         endif
1816       endif
1817       call chainbuild_cart
1818       call kinetic(EK)
1819       if (tbf) then
1820         call verlet_bath
1821       endif
1822       kinetic_T=2.0d0/(dimen3*Rb)*EK
1823       if(me.eq.king.or..not.out1file)then
1824        call cartprint
1825        call intout
1826       endif
1827 #ifdef MPI
1828       tt0=MPI_Wtime()
1829 #else
1830       tt0=tcpu()
1831 #endif
1832       call zerograd
1833       call etotal(potEcomp)
1834       if (large) call enerprint(potEcomp)
1835 #ifdef TIMING_ENE
1836 #ifdef MPI
1837       t_etotal=t_etotal+MPI_Wtime()-tt0
1838 #else
1839       t_etotal=t_etotal+tcpu()-tt0
1840 #endif
1841 #endif
1842       potE=potEcomp(0)
1843       call cartgrad
1844       call lagrangian
1845       call max_accel
1846       if (amax*d_time .gt. dvmax) then
1847         d_time=d_time*dvmax/amax
1848         if(me.eq.king.or..not.out1file) write (iout,*)
1849      &   "Time step reduced to",d_time,
1850      &   " because of too large initial acceleration."
1851       endif
1852       if(me.eq.king.or..not.out1file)then
1853        write(iout,*) "Potential energy and its components"
1854        call enerprint(potEcomp)
1855        write(iout,*) (potEcomp(i),i=0,n_ene)
1856       endif
1857       potE=potEcomp(0)-potEcomp(20)
1858       totE=EK+potE
1859       itime=0
1860       if (ntwe.ne.0) call statout(itime)
1861       if(me.eq.king.or..not.out1file)
1862      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1863      &   " Kinetic energy",EK," Potential energy",potE,
1864      &   " Total energy",totE," Maximum acceleration ",
1865      &   amax
1866       if (large) then
1867         write (iout,*) "Initial coordinates"
1868         do i=1,nres
1869           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1870      &    (c(j,i+nres),j=1,3)
1871         enddo
1872         write (iout,*) "Initial dC"
1873         do i=0,nres
1874           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1875      &    (dc(j,i+nres),j=1,3)
1876         enddo
1877         write (iout,*) "Initial velocities"
1878         write (iout,"(13x,' backbone ',23x,' side chain')")
1879         do i=0,nres
1880           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1881      &    (d_t(j,i+nres),j=1,3)
1882         enddo
1883         write (iout,*) "Initial accelerations"
1884         do i=0,nres
1885 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1886           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1887      &    (d_a(j,i+nres),j=1,3)
1888         enddo
1889       endif
1890       do i=0,2*nres
1891         do j=1,3
1892           dc_old(j,i)=dc(j,i)
1893           d_t_old(j,i)=d_t(j,i)
1894           d_a_old(j,i)=d_a(j,i)
1895         enddo
1896 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1897       enddo
1898       if (RESPA) then
1899 #ifdef MPI
1900         tt0 =MPI_Wtime()
1901 #else
1902         tt0 = tcpu()
1903 #endif
1904         call zerograd
1905         call etotal_short(energia_short)
1906         if (large) call enerprint(potEcomp)
1907 #ifdef TIMING_ENE
1908 #ifdef MPI
1909         t_eshort=t_eshort+MPI_Wtime()-tt0
1910 #else
1911         t_eshort=t_eshort+tcpu()-tt0
1912 #endif
1913 #endif
1914         call cartgrad
1915         call lagrangian
1916         if(.not.out1file .and. large) then
1917           write (iout,*) "energia_long",energia_long(0),
1918      &     " energia_short",energia_short(0),
1919      &     " total",energia_long(0)+energia_short(0)
1920           write (iout,*) "Initial fast-force accelerations"
1921           do i=0,nres
1922             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1923      &      (d_a(j,i+nres),j=1,3)
1924           enddo
1925         endif
1926 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1927         do i=0,2*nres
1928           do j=1,3
1929             d_a_short(j,i)=d_a(j,i)
1930           enddo
1931         enddo
1932 #ifdef MPI
1933         tt0=MPI_Wtime()
1934 #else
1935         tt0=tcpu()
1936 #endif
1937         call zerograd
1938         call etotal_long(energia_long)
1939         if (large) call enerprint(potEcomp)
1940 #ifdef TIMING_ENE
1941 #ifdef MPI
1942         t_elong=t_elong+MPI_Wtime()-tt0
1943 #else
1944         t_elong=t_elong+tcpu()-tt0
1945 #endif
1946 #endif
1947         call cartgrad
1948         call lagrangian
1949         if(.not.out1file .and. large) then
1950           write (iout,*) "energia_long",energia_long(0)
1951           write (iout,*) "Initial slow-force accelerations"
1952           do i=0,nres
1953             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1954      &      (d_a(j,i+nres),j=1,3)
1955           enddo
1956         endif
1957 #ifdef MPI
1958         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1959 #else
1960         t_enegrad=t_enegrad+tcpu()-tt0
1961 #endif
1962       endif
1963       return
1964       end
1965 c-----------------------------------------------------------
1966       subroutine random_vel
1967       implicit real*8 (a-h,o-z)
1968       include 'DIMENSIONS'
1969       include 'COMMON.CONTROL'
1970       include 'COMMON.VAR'
1971       include 'COMMON.MD'
1972       include 'COMMON.LAGRANGE'
1973 #ifndef LANG0
1974       include 'COMMON.LANGEVIN'
1975 #else
1976       include 'COMMON.LANGEVIN.lang0'
1977 #endif
1978       include 'COMMON.CHAIN'
1979       include 'COMMON.DERIV'
1980       include 'COMMON.GEO'
1981       include 'COMMON.LOCAL'
1982       include 'COMMON.INTERACT'
1983       include 'COMMON.IOUNITS'
1984       include 'COMMON.NAMES'
1985       include 'COMMON.TIME1'
1986       double precision xv,sigv,lowb,highb,vec_afm(3)
1987 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1988 c First generate velocities in the eigenspace of the G matrix
1989 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1990 c      call flush(iout)
1991       xv=0.0d0
1992       ii=0
1993       do i=1,dimen
1994         do k=1,3
1995           ii=ii+1
1996           sigv=dsqrt((Rb*t_bath)/geigen(i))
1997           lowb=-5*sigv
1998           highb=5*sigv
1999           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2000
2001 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2002 c     &      " d_t_work_new",d_t_work_new(ii)
2003         enddo
2004       enddo
2005 C       if (SELFGUIDE.gt.0) then
2006 C       distance=0.0
2007 C       do j=1,3
2008 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2009 C       distance=distance+vec_afm(j)**2
2010 C       enddo
2011 C       distance=dsqrt(distance)
2012 C       do j=1,3
2013 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2014 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2015 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2016 C     &    d_t_work_new(j+(afmend-1)*3)
2017 C       enddo
2018
2019 C       endif
2020
2021 c diagnostics
2022 c      Ek1=0.0d0
2023 c      ii=0
2024 c      do i=1,dimen
2025 c        do k=1,3
2026 c          ii=ii+1
2027 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2028 c        enddo
2029 c      enddo
2030 c      write (iout,*) "Ek from eigenvectors",Ek1
2031 c end diagnostics
2032 c Transform velocities to UNRES coordinate space
2033       do k=0,2
2034         do i=1,dimen
2035           ind=(i-1)*3+k+1
2036           d_t_work(ind)=0.0d0
2037           do j=1,dimen
2038             d_t_work(ind)=d_t_work(ind)
2039      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2040           enddo
2041 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2042 c          call flush(iout)
2043         enddo
2044       enddo
2045 c Transfer to the d_t vector
2046       do j=1,3
2047         d_t(j,0)=d_t_work(j)
2048       enddo
2049       ind=3
2050       do i=nnt,nct-1
2051         do j=1,3
2052           ind=ind+1
2053           d_t(j,i)=d_t_work(ind)
2054         enddo
2055       enddo
2056       do i=nnt,nct
2057         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2058           do j=1,3
2059             ind=ind+1
2060             d_t(j,i+nres)=d_t_work(ind)
2061           enddo
2062         endif
2063       enddo
2064 c      call kinetic(EK)
2065 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2066 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2067 c      call flush(iout)
2068       return
2069       end
2070 #ifndef LANG0
2071 c-----------------------------------------------------------
2072       subroutine sd_verlet_p_setup
2073 c Sets up the parameters of stochastic Verlet algorithm
2074       implicit real*8 (a-h,o-z)
2075       include 'DIMENSIONS'
2076 #ifdef MPI
2077       include 'mpif.h'
2078 #endif
2079       include 'COMMON.CONTROL'
2080       include 'COMMON.VAR'
2081       include 'COMMON.MD'
2082 #ifndef LANG0
2083       include 'COMMON.LANGEVIN'
2084 #else
2085       include 'COMMON.LANGEVIN.lang0'
2086 #endif
2087       include 'COMMON.CHAIN'
2088       include 'COMMON.DERIV'
2089       include 'COMMON.GEO'
2090       include 'COMMON.LOCAL'
2091       include 'COMMON.INTERACT'
2092       include 'COMMON.IOUNITS'
2093       include 'COMMON.NAMES'
2094       include 'COMMON.TIME1'
2095       double precision emgdt(MAXRES6),
2096      & pterm,vterm,rho,rhoc,vsig,
2097      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2098      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2099      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2100       logical lprn /.false./
2101       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2102       double precision ktm
2103 #ifdef MPI
2104       tt0 = MPI_Wtime()
2105 #else
2106       tt0 = tcpu()
2107 #endif
2108 c
2109 c AL 8/17/04 Code adapted from tinker
2110 c
2111 c Get the frictional and random terms for stochastic dynamics in the
2112 c eigenspace of mass-scaled UNRES friction matrix
2113 c
2114       do i = 1, dimen
2115             gdt = fricgam(i) * d_time
2116 c
2117 c Stochastic dynamics reduces to simple MD for zero friction
2118 c
2119             if (gdt .le. zero) then
2120                pfric_vec(i) = 1.0d0
2121                vfric_vec(i) = d_time
2122                afric_vec(i) = 0.5d0 * d_time * d_time
2123                prand_vec(i) = 0.0d0
2124                vrand_vec1(i) = 0.0d0
2125                vrand_vec2(i) = 0.0d0
2126 c
2127 c Analytical expressions when friction coefficient is large
2128 c
2129             else
2130                if (gdt .ge. gdt_radius) then
2131                   egdt = dexp(-gdt)
2132                   pfric_vec(i) = egdt
2133                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2134                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2135                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2136                   vterm = 1.0d0 - egdt**2
2137                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2138 c
2139 c Use series expansions when friction coefficient is small
2140 c
2141                else
2142                   gdt2 = gdt * gdt
2143                   gdt3 = gdt * gdt2
2144                   gdt4 = gdt2 * gdt2
2145                   gdt5 = gdt2 * gdt3
2146                   gdt6 = gdt3 * gdt3
2147                   gdt7 = gdt3 * gdt4
2148                   gdt8 = gdt4 * gdt4
2149                   gdt9 = gdt4 * gdt5
2150                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2151      &                          - gdt5/120.0d0 + gdt6/720.0d0
2152      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2153      &                          - gdt9/362880.0d0) / fricgam(i)**2
2154                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2155                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2156                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2157      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2158      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2159      &                       + 127.0d0*gdt9/90720.0d0
2160                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2161      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2162      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2163      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2164                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2165      &                       - 17.0d0*gdt2/1280.0d0
2166      &                       + 17.0d0*gdt3/6144.0d0
2167      &                       + 40967.0d0*gdt4/34406400.0d0
2168      &                       - 57203.0d0*gdt5/275251200.0d0
2169      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2170                end if
2171 c
2172 c Compute the scaling factors of random terms for the nonzero friction case
2173 c
2174                ktm = 0.5d0*d_time/fricgam(i)
2175                psig = dsqrt(ktm*pterm) / fricgam(i)
2176                vsig = dsqrt(ktm*vterm)
2177                rhoc = dsqrt(1.0d0 - rho*rho)
2178                prand_vec(i) = psig
2179                vrand_vec1(i) = vsig * rho
2180                vrand_vec2(i) = vsig * rhoc
2181             end if
2182       end do
2183       if (lprn) then
2184       write (iout,*)
2185      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2186      &  " vrand_vec2"
2187       do i=1,dimen
2188         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2189      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2190       enddo
2191       endif
2192 c
2193 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2194 c
2195 #ifndef   LANG0
2196       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2197       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2198       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2199       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2200       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2201       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2202 #endif
2203 #ifdef MPI
2204       t_sdsetup=t_sdsetup+MPI_Wtime()
2205 #else
2206       t_sdsetup=t_sdsetup+tcpu()-tt0
2207 #endif
2208       return
2209       end
2210 c-------------------------------------------------------------
2211       subroutine eigtransf1(n,ndim,ab,d,c)
2212       implicit none
2213       integer n,ndim
2214       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2215       integer i,j,k
2216       do i=1,n
2217         do j=1,n
2218           c(i,j)=0.0d0
2219           do k=1,n
2220             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2221           enddo
2222         enddo
2223       enddo
2224       return
2225       end
2226 c-------------------------------------------------------------
2227       subroutine eigtransf(n,ndim,a,b,d,c)
2228       implicit none
2229       integer n,ndim
2230       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2231       integer i,j,k
2232       do i=1,n
2233         do j=1,n
2234           c(i,j)=0.0d0
2235           do k=1,n
2236             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2237           enddo
2238         enddo
2239       enddo
2240       return
2241       end
2242 c-------------------------------------------------------------
2243       subroutine sd_verlet1
2244 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2245       implicit real*8 (a-h,o-z)
2246       include 'DIMENSIONS'
2247       include 'COMMON.CONTROL'
2248       include 'COMMON.VAR'
2249       include 'COMMON.MD'
2250       include 'COMMON.LAGRANGE'
2251 #ifndef LANG0
2252       include 'COMMON.LANGEVIN'
2253 #else
2254       include 'COMMON.LANGEVIN.lang0'
2255 #endif
2256       include 'COMMON.CHAIN'
2257       include 'COMMON.DERIV'
2258       include 'COMMON.GEO'
2259       include 'COMMON.LOCAL'
2260       include 'COMMON.INTERACT'
2261       include 'COMMON.IOUNITS'
2262       include 'COMMON.NAMES'
2263       double precision stochforcvec(MAXRES6)
2264       common /stochcalc/ stochforcvec
2265       logical lprn /.false./
2266
2267 c      write (iout,*) "dc_old"
2268 c      do i=0,nres
2269 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2270 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2271 c      enddo
2272       do j=1,3
2273         dc_work(j)=dc_old(j,0)
2274         d_t_work(j)=d_t_old(j,0)
2275         d_a_work(j)=d_a_old(j,0)
2276       enddo
2277       ind=3
2278       do i=nnt,nct-1
2279         do j=1,3
2280           dc_work(ind+j)=dc_old(j,i)
2281           d_t_work(ind+j)=d_t_old(j,i)
2282           d_a_work(ind+j)=d_a_old(j,i)
2283         enddo
2284         ind=ind+3
2285       enddo
2286       do i=nnt,nct
2287         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2288           do j=1,3
2289             dc_work(ind+j)=dc_old(j,i+nres)
2290             d_t_work(ind+j)=d_t_old(j,i+nres)
2291             d_a_work(ind+j)=d_a_old(j,i+nres)
2292           enddo
2293           ind=ind+3
2294         endif
2295       enddo
2296 #ifndef LANG0
2297       if (lprn) then
2298       write (iout,*)
2299      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2300      &  " vrand_mat2"
2301       do i=1,dimen
2302         do j=1,dimen
2303           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2304      &      vfric_mat(i,j),afric_mat(i,j),
2305      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2306         enddo
2307       enddo
2308       endif
2309       do i=1,dimen
2310         ddt1=0.0d0
2311         ddt2=0.0d0
2312         do j=1,dimen
2313           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2314      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2315           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2316           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2317         enddo
2318         d_t_work_new(i)=ddt1+0.5d0*ddt2
2319         d_t_work(i)=ddt1+ddt2
2320       enddo
2321 #endif
2322       do j=1,3
2323         dc(j,0)=dc_work(j)
2324         d_t(j,0)=d_t_work(j)
2325       enddo
2326       ind=3
2327       do i=nnt,nct-1
2328         do j=1,3
2329           dc(j,i)=dc_work(ind+j)
2330           d_t(j,i)=d_t_work(ind+j)
2331         enddo
2332         ind=ind+3
2333       enddo
2334       do i=nnt,nct
2335         if (itype(i).ne.10) then
2336           inres=i+nres
2337           do j=1,3
2338             dc(j,inres)=dc_work(ind+j)
2339             d_t(j,inres)=d_t_work(ind+j)
2340           enddo
2341           ind=ind+3
2342         endif
2343       enddo
2344       return
2345       end
2346 c--------------------------------------------------------------------------
2347       subroutine sd_verlet2
2348 c  Calculating the adjusted velocities for accelerations
2349       implicit real*8 (a-h,o-z)
2350       include 'DIMENSIONS'
2351       include 'COMMON.CONTROL'
2352       include 'COMMON.VAR'
2353       include 'COMMON.MD'
2354       include 'COMMON.LAGRANGE'
2355 #ifndef LANG0
2356       include 'COMMON.LANGEVIN'
2357 #else
2358       include 'COMMON.LANGEVIN.lang0'
2359 #endif
2360       include 'COMMON.CHAIN'
2361       include 'COMMON.DERIV'
2362       include 'COMMON.GEO'
2363       include 'COMMON.LOCAL'
2364       include 'COMMON.INTERACT'
2365       include 'COMMON.IOUNITS'
2366       include 'COMMON.NAMES'
2367       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2368       common /stochcalc/ stochforcvec
2369 c
2370 c Compute the stochastic forces which contribute to velocity change
2371 c
2372       call stochastic_force(stochforcvecV)
2373
2374 #ifndef LANG0
2375       do i=1,dimen
2376         ddt1=0.0d0
2377         ddt2=0.0d0
2378         do j=1,dimen
2379           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2380           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2381      &     vrand_mat2(i,j)*stochforcvecV(j)
2382         enddo
2383         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2384       enddo
2385 #endif
2386       do j=1,3
2387         d_t(j,0)=d_t_work(j)
2388       enddo
2389       ind=3
2390       do i=nnt,nct-1
2391         do j=1,3
2392           d_t(j,i)=d_t_work(ind+j)
2393         enddo
2394         ind=ind+3
2395       enddo
2396       do i=nnt,nct
2397         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2398           inres=i+nres
2399           do j=1,3
2400             d_t(j,inres)=d_t_work(ind+j)
2401           enddo
2402           ind=ind+3
2403         endif
2404       enddo
2405       return
2406       end
2407 c-----------------------------------------------------------
2408       subroutine sd_verlet_ciccotti_setup
2409 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2410 c version
2411       implicit real*8 (a-h,o-z)
2412       include 'DIMENSIONS'
2413 #ifdef MPI
2414       include 'mpif.h'
2415 #endif
2416       include 'COMMON.CONTROL'
2417       include 'COMMON.VAR'
2418       include 'COMMON.MD'
2419       include 'COMMON.LAGRANGE'
2420 #ifndef LANG0
2421       include 'COMMON.LANGEVIN'
2422 #else
2423       include 'COMMON.LANGEVIN.lang0'
2424 #endif
2425       include 'COMMON.CHAIN'
2426       include 'COMMON.DERIV'
2427       include 'COMMON.GEO'
2428       include 'COMMON.LOCAL'
2429       include 'COMMON.INTERACT'
2430       include 'COMMON.IOUNITS'
2431       include 'COMMON.NAMES'
2432       include 'COMMON.TIME1'
2433       double precision emgdt(MAXRES6),
2434      & pterm,vterm,rho,rhoc,vsig,
2435      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2436      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2437      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2438       logical lprn /.false./
2439       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2440       double precision ktm
2441 #ifdef MPI
2442       tt0 = MPI_Wtime()
2443 #else
2444       tt0 = tcpu()
2445 #endif
2446 c
2447 c AL 8/17/04 Code adapted from tinker
2448 c
2449 c Get the frictional and random terms for stochastic dynamics in the
2450 c eigenspace of mass-scaled UNRES friction matrix
2451 c
2452       do i = 1, dimen
2453             write (iout,*) "i",i," fricgam",fricgam(i)
2454             gdt = fricgam(i) * d_time
2455 c
2456 c Stochastic dynamics reduces to simple MD for zero friction
2457 c
2458             if (gdt .le. zero) then
2459                pfric_vec(i) = 1.0d0
2460                vfric_vec(i) = d_time
2461                afric_vec(i) = 0.5d0*d_time*d_time
2462                prand_vec(i) = afric_vec(i)
2463                vrand_vec2(i) = vfric_vec(i)
2464 c
2465 c Analytical expressions when friction coefficient is large
2466 c
2467             else
2468                egdt = dexp(-gdt)
2469                pfric_vec(i) = egdt
2470                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2471                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2472                prand_vec(i) = afric_vec(i)
2473                vrand_vec2(i) = vfric_vec(i)
2474 c
2475 c Compute the scaling factors of random terms for the nonzero friction case
2476 c
2477 c               ktm = 0.5d0*d_time/fricgam(i)
2478 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2479 c               vsig = dsqrt(ktm*vterm)
2480 c               prand_vec(i) = psig*afric_vec(i)
2481 c               vrand_vec2(i) = vsig*vfric_vec(i)
2482             end if
2483       end do
2484       if (lprn) then
2485       write (iout,*)
2486      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2487      &  " vrand_vec2"
2488       do i=1,dimen
2489         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2490      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2491       enddo
2492       endif
2493 c
2494 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2495 c
2496       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2497       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2498       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2499       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2500       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2501 #ifdef MPI
2502       t_sdsetup=t_sdsetup+MPI_Wtime()
2503 #else
2504       t_sdsetup=t_sdsetup+tcpu()-tt0
2505 #endif
2506       return
2507       end
2508 c-------------------------------------------------------------
2509       subroutine sd_verlet1_ciccotti
2510 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2511       implicit real*8 (a-h,o-z)
2512       include 'DIMENSIONS'
2513 #ifdef MPI
2514       include 'mpif.h'
2515 #endif
2516       include 'COMMON.CONTROL'
2517       include 'COMMON.VAR'
2518       include 'COMMON.MD'
2519       include 'COMMON.LAGRANGE'
2520 #ifndef LANG0
2521       include 'COMMON.LANGEVIN'
2522 #else
2523       include 'COMMON.LANGEVIN.lang0'
2524 #endif
2525       include 'COMMON.CHAIN'
2526       include 'COMMON.DERIV'
2527       include 'COMMON.GEO'
2528       include 'COMMON.LOCAL'
2529       include 'COMMON.INTERACT'
2530       include 'COMMON.IOUNITS'
2531       include 'COMMON.NAMES'
2532       double precision stochforcvec(MAXRES6)
2533       common /stochcalc/ stochforcvec
2534       logical lprn /.false./
2535
2536 c      write (iout,*) "dc_old"
2537 c      do i=0,nres
2538 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2539 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2540 c      enddo
2541       do j=1,3
2542         dc_work(j)=dc_old(j,0)
2543         d_t_work(j)=d_t_old(j,0)
2544         d_a_work(j)=d_a_old(j,0)
2545       enddo
2546       ind=3
2547       do i=nnt,nct-1
2548         do j=1,3
2549           dc_work(ind+j)=dc_old(j,i)
2550           d_t_work(ind+j)=d_t_old(j,i)
2551           d_a_work(ind+j)=d_a_old(j,i)
2552         enddo
2553         ind=ind+3
2554       enddo
2555       do i=nnt,nct
2556         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2557           do j=1,3
2558             dc_work(ind+j)=dc_old(j,i+nres)
2559             d_t_work(ind+j)=d_t_old(j,i+nres)
2560             d_a_work(ind+j)=d_a_old(j,i+nres)
2561           enddo
2562           ind=ind+3
2563         endif
2564       enddo
2565
2566 #ifndef LANG0
2567       if (lprn) then
2568       write (iout,*)
2569      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2570      &  " vrand_mat2"
2571       do i=1,dimen
2572         do j=1,dimen
2573           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2574      &      vfric_mat(i,j),afric_mat(i,j),
2575      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2576         enddo
2577       enddo
2578       endif
2579       do i=1,dimen
2580         ddt1=0.0d0
2581         ddt2=0.0d0
2582         do j=1,dimen
2583           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2584      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2585           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2586           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2587         enddo
2588         d_t_work_new(i)=ddt1+0.5d0*ddt2
2589         d_t_work(i)=ddt1+ddt2
2590       enddo
2591 #endif
2592       do j=1,3
2593         dc(j,0)=dc_work(j)
2594         d_t(j,0)=d_t_work(j)
2595       enddo
2596       ind=3
2597       do i=nnt,nct-1
2598         do j=1,3
2599           dc(j,i)=dc_work(ind+j)
2600           d_t(j,i)=d_t_work(ind+j)
2601         enddo
2602         ind=ind+3
2603       enddo
2604       do i=nnt,nct
2605         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2606           inres=i+nres
2607           do j=1,3
2608             dc(j,inres)=dc_work(ind+j)
2609             d_t(j,inres)=d_t_work(ind+j)
2610           enddo
2611           ind=ind+3
2612         endif
2613       enddo
2614       return
2615       end
2616 c--------------------------------------------------------------------------
2617       subroutine sd_verlet2_ciccotti
2618 c  Calculating the adjusted velocities for accelerations
2619       implicit real*8 (a-h,o-z)
2620       include 'DIMENSIONS'
2621       include 'COMMON.CONTROL'
2622       include 'COMMON.VAR'
2623       include 'COMMON.MD'
2624       iclude 'COMMON.LAGRANGE'
2625 #ifndef LANG0
2626       include 'COMMON.LANGEVIN'
2627 #else
2628       include 'COMMON.LANGEVIN.lang0'
2629 #endif
2630       include 'COMMON.CHAIN'
2631       include 'COMMON.DERIV'
2632       include 'COMMON.GEO'
2633       include 'COMMON.LOCAL'
2634       include 'COMMON.INTERACT'
2635       include 'COMMON.IOUNITS'
2636       include 'COMMON.NAMES'
2637       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2638       common /stochcalc/ stochforcvec
2639 c
2640 c Compute the stochastic forces which contribute to velocity change
2641 c
2642       call stochastic_force(stochforcvecV)
2643 #ifndef LANG0
2644       do i=1,dimen
2645         ddt1=0.0d0
2646         ddt2=0.0d0
2647         do j=1,dimen
2648
2649           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2650 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2651           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2652         enddo
2653         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2654       enddo
2655 #endif
2656       do j=1,3
2657         d_t(j,0)=d_t_work(j)
2658       enddo
2659       ind=3
2660       do i=nnt,nct-1
2661         do j=1,3
2662           d_t(j,i)=d_t_work(ind+j)
2663         enddo
2664         ind=ind+3
2665       enddo
2666       do i=nnt,nct
2667         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2668           inres=i+nres
2669           do j=1,3
2670             d_t(j,inres)=d_t_work(ind+j)
2671           enddo
2672           ind=ind+3
2673         endif
2674       enddo
2675       return
2676       end
2677 #endif