source/unres/src_MD-M/MD_A-MTS.F

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