source/unres/src_MD-M-SAXS.safe/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           if (itype(i).ne.ntyp1) 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       endif
1653       call inertia_tensor
1654 c  Getting the potential energy and forces and velocities and accelerations
1655       call vcm_vel(vcm)
1656       write (iout,*) "velocity of the center of the mass:"
1657       write (iout,*) (vcm(j),j=1,3)
1658       call flush(iout)
1659       do j=1,3
1660         d_t(j,0)=d_t(j,0)-vcm(j)
1661       enddo
1662 c Removing the velocity of the center of mass
1663       call vcm_vel(vcm)
1664       if(me.eq.king.or..not.out1file)then
1665        write (iout,*) "vcm right after adjustment:"
1666        write (iout,*) (vcm(j),j=1,3)
1667       endif
1668       call flush(iout)
1669       if (.not.rest) then
1670          if(iranconf.ne.0) then
1671 c 8/22/17 AL Loop to produce a low-energy random conformation
1672           do iranmin=1,10
1673           call chainbuild
1674           if (overlapsc) then
1675            print *, 'Calling OVERLAP_SC'
1676            call overlap_sc(fail)
1677           endif
1678
1679           if (searchsc) then
1680            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1681            print *,'SC_move',nft_sc,etot
1682            if(me.eq.king.or..not.out1file)
1683      &      write(iout,*) 'SC_move',nft_sc,etot
1684           endif
1685
1686           if(dccart)then
1687            print *, 'Calling MINIM_DC'
1688            call minim_dc(etot,iretcode,nfun)
1689           else
1690            call geom_to_var(nvar,varia)
1691            print *,'Calling MINIMIZE.'
1692            call minimize(etot,varia,iretcode,nfun)
1693            call var_to_geom(nvar,varia)
1694           endif
1695           if(me.eq.king.or..not.out1file)
1696      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1697           if (isnan(etot) .or. etot.gt.1.0d4) then
1698             write (iout,*) "Energy too large",etot,
1699      &        " trying another random conformation"
1700             do itrial=1,100
1701               itmp=1
1702               call gen_rand_conf(itmp,*30)
1703               goto 40
1704    30         write (iout,*) 'Failed to generate random conformation',
1705      &          ', itrial=',itrial
1706               write (*,*) 'Processor:',me,
1707      &          ' Failed to generate random conformation',
1708      &          ' itrial=',itrial
1709               call intout
1710
1711 #ifdef AIX
1712               call flush_(iout)
1713 #else
1714               call flush(iout)
1715 #endif
1716             enddo
1717             write (iout,'(a,i3,a)') 'Processor:',me,
1718      &        ' error in generating random conformation.'
1719             write (*,'(a,i3,a)') 'Processor:',me,
1720      &        ' error in generating random conformation.'
1721             call flush(iout)
1722 #ifdef MPI
1723             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1724 #else
1725             stop
1726 #endif
1727    40       continue
1728           else
1729             goto 44
1730           endif
1731           enddo
1732           write (iout,'(a,i3,a)') 'Processor:',me,
1733      &        ' failed to generate a low-energy random conformation.'
1734             write (*,'(a,i3,a)') 'Processor:',me,
1735      &        ' failed to generate a low-energy random conformation.'
1736             call flush(iout)
1737 #ifdef MPI
1738             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1739 #else
1740             stop
1741 #endif
1742    44     continue
1743          else if (indpdb.gt.0) then
1744 C 8/22/17 AL Minimize initial PDB structure
1745           if(dccart)then
1746            print *, 'Calling MINIM_DC'
1747            call minim_dc(etot,iretcode,nfun)
1748           else
1749            call geom_to_var(nvar,varia)
1750            print *,'Calling MINIMIZE.'
1751            call minimize(etot,varia,iretcode,nfun)
1752            call var_to_geom(nvar,varia)
1753           endif
1754           if(me.eq.king.or..not.out1file)
1755      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1756          endif
1757       endif
1758       call chainbuild_cart
1759       call kinetic(EK)
1760       if (tbf) then
1761         call verlet_bath
1762       endif
1763       kinetic_T=2.0d0/(dimen3*Rb)*EK
1764       if(me.eq.king.or..not.out1file)then
1765        call cartprint
1766        call intout
1767       endif
1768 #ifdef MPI
1769       tt0=MPI_Wtime()
1770 #else
1771       tt0=tcpu()
1772 #endif
1773       call zerograd
1774       call etotal(potEcomp)
1775       if (large) call enerprint(potEcomp)
1776 #ifdef TIMING_ENE
1777 #ifdef MPI
1778       t_etotal=t_etotal+MPI_Wtime()-tt0
1779 #else
1780       t_etotal=t_etotal+tcpu()-tt0
1781 #endif
1782 #endif
1783       potE=potEcomp(0)
1784       call cartgrad
1785       call lagrangian
1786       call max_accel
1787       if (amax*d_time .gt. dvmax) then
1788         d_time=d_time*dvmax/amax
1789         if(me.eq.king.or..not.out1file) write (iout,*)
1790      &   "Time step reduced to",d_time,
1791      &   " because of too large initial acceleration."
1792       endif
1793       if(me.eq.king.or..not.out1file)then
1794        write(iout,*) "Potential energy and its components"
1795        call enerprint(potEcomp)
1796        write(iout,*) (potEcomp(i),i=0,n_ene)
1797       endif
1798       potE=potEcomp(0)-potEcomp(20)
1799       totE=EK+potE
1800       itime=0
1801       if (ntwe.ne.0) call statout(itime)
1802       if(me.eq.king.or..not.out1file)
1803      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1804      &   " Kinetic energy",EK," Potential energy",potE,
1805      &   " Total energy",totE," Maximum acceleration ",
1806      &   amax
1807       if (large) then
1808         write (iout,*) "Initial coordinates"
1809         do i=1,nres
1810           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1811      &    (c(j,i+nres),j=1,3)
1812         enddo
1813         write (iout,*) "Initial dC"
1814         do i=0,nres
1815           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1816      &    (dc(j,i+nres),j=1,3)
1817         enddo
1818         write (iout,*) "Initial velocities"
1819         write (iout,"(13x,' backbone ',23x,' side chain')")
1820         do i=0,nres
1821           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1822      &    (d_t(j,i+nres),j=1,3)
1823         enddo
1824         write (iout,*) "Initial accelerations"
1825         do i=0,nres
1826 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1827           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1828      &    (d_a(j,i+nres),j=1,3)
1829         enddo
1830       endif
1831       do i=0,2*nres
1832         do j=1,3
1833           dc_old(j,i)=dc(j,i)
1834           d_t_old(j,i)=d_t(j,i)
1835           d_a_old(j,i)=d_a(j,i)
1836         enddo
1837 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1838       enddo
1839       if (RESPA) then
1840 #ifdef MPI
1841         tt0 =MPI_Wtime()
1842 #else
1843         tt0 = tcpu()
1844 #endif
1845         call zerograd
1846         call etotal_short(energia_short)
1847         if (large) call enerprint(potEcomp)
1848 #ifdef TIMING_ENE
1849 #ifdef MPI
1850         t_eshort=t_eshort+MPI_Wtime()-tt0
1851 #else
1852         t_eshort=t_eshort+tcpu()-tt0
1853 #endif
1854 #endif
1855         call cartgrad
1856         call lagrangian
1857         if(.not.out1file .and. large) then
1858           write (iout,*) "energia_long",energia_long(0),
1859      &     " energia_short",energia_short(0),
1860      &     " total",energia_long(0)+energia_short(0)
1861           write (iout,*) "Initial fast-force accelerations"
1862           do i=0,nres
1863             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1864      &      (d_a(j,i+nres),j=1,3)
1865           enddo
1866         endif
1867 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1868         do i=0,2*nres
1869           do j=1,3
1870             d_a_short(j,i)=d_a(j,i)
1871           enddo
1872         enddo
1873 #ifdef MPI
1874         tt0=MPI_Wtime()
1875 #else
1876         tt0=tcpu()
1877 #endif
1878         call zerograd
1879         call etotal_long(energia_long)
1880         if (large) call enerprint(potEcomp)
1881 #ifdef TIMING_ENE
1882 #ifdef MPI
1883         t_elong=t_elong+MPI_Wtime()-tt0
1884 #else
1885         t_elong=t_elong+tcpu()-tt0
1886 #endif
1887 #endif
1888         call cartgrad
1889         call lagrangian
1890         if(.not.out1file .and. large) then
1891           write (iout,*) "energia_long",energia_long(0)
1892           write (iout,*) "Initial slow-force accelerations"
1893           do i=0,nres
1894             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1895      &      (d_a(j,i+nres),j=1,3)
1896           enddo
1897         endif
1898 #ifdef MPI
1899         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1900 #else
1901         t_enegrad=t_enegrad+tcpu()-tt0
1902 #endif
1903       endif
1904       return
1905       end
1906 c-----------------------------------------------------------
1907       subroutine random_vel
1908       implicit real*8 (a-h,o-z)
1909       include 'DIMENSIONS'
1910       include 'COMMON.CONTROL'
1911       include 'COMMON.VAR'
1912       include 'COMMON.MD'
1913 #ifndef LANG0
1914       include 'COMMON.LANGEVIN'
1915 #else
1916       include 'COMMON.LANGEVIN.lang0'
1917 #endif
1918       include 'COMMON.CHAIN'
1919       include 'COMMON.DERIV'
1920       include 'COMMON.GEO'
1921       include 'COMMON.LOCAL'
1922       include 'COMMON.INTERACT'
1923       include 'COMMON.IOUNITS'
1924       include 'COMMON.NAMES'
1925       include 'COMMON.TIME1'
1926       double precision xv,sigv,lowb,highb,vec_afm(3)
1927 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1928 c First generate velocities in the eigenspace of the G matrix
1929 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1930 c      call flush(iout)
1931       xv=0.0d0
1932       ii=0
1933       do i=1,dimen
1934         do k=1,3
1935           ii=ii+1
1936           sigv=dsqrt((Rb*t_bath)/geigen(i))
1937           lowb=-5*sigv
1938           highb=5*sigv
1939           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1940
1941 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1942 c     &      " d_t_work_new",d_t_work_new(ii)
1943         enddo
1944       enddo
1945 C       if (SELFGUIDE.gt.0) then
1946 C       distance=0.0
1947 C       do j=1,3
1948 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1949 C       distance=distance+vec_afm(j)**2
1950 C       enddo
1951 C       distance=dsqrt(distance)
1952 C       do j=1,3
1953 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1954 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1955 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1956 C     &    d_t_work_new(j+(afmend-1)*3)
1957 C       enddo
1958
1959 C       endif
1960
1961 c diagnostics
1962 c      Ek1=0.0d0
1963 c      ii=0
1964 c      do i=1,dimen
1965 c        do k=1,3
1966 c          ii=ii+1
1967 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1968 c        enddo
1969 c      enddo
1970 c      write (iout,*) "Ek from eigenvectors",Ek1
1971 c end diagnostics
1972 c Transform velocities to UNRES coordinate space
1973       do k=0,2
1974         do i=1,dimen
1975           ind=(i-1)*3+k+1
1976           d_t_work(ind)=0.0d0
1977           do j=1,dimen
1978             d_t_work(ind)=d_t_work(ind)
1979      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1980           enddo
1981 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1982 c          call flush(iout)
1983         enddo
1984       enddo
1985 c Transfer to the d_t vector
1986       do j=1,3
1987         d_t(j,0)=d_t_work(j)
1988       enddo
1989       ind=3
1990       do i=nnt,nct-1
1991         do j=1,3
1992           ind=ind+1
1993           d_t(j,i)=d_t_work(ind)
1994         enddo
1995       enddo
1996       do i=nnt,nct
1997         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1998           do j=1,3
1999             ind=ind+1
2000             d_t(j,i+nres)=d_t_work(ind)
2001           enddo
2002         endif
2003       enddo
2004 c      call kinetic(EK)
2005 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2006 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2007 c      call flush(iout)
2008       return
2009       end
2010 #ifndef LANG0
2011 c-----------------------------------------------------------
2012       subroutine sd_verlet_p_setup
2013 c Sets up the parameters of stochastic Verlet algorithm
2014       implicit real*8 (a-h,o-z)
2015       include 'DIMENSIONS'
2016 #ifdef MPI
2017       include 'mpif.h'
2018 #endif
2019       include 'COMMON.CONTROL'
2020       include 'COMMON.VAR'
2021       include 'COMMON.MD'
2022 #ifndef LANG0
2023       include 'COMMON.LANGEVIN'
2024 #else
2025       include 'COMMON.LANGEVIN.lang0'
2026 #endif
2027       include 'COMMON.CHAIN'
2028       include 'COMMON.DERIV'
2029       include 'COMMON.GEO'
2030       include 'COMMON.LOCAL'
2031       include 'COMMON.INTERACT'
2032       include 'COMMON.IOUNITS'
2033       include 'COMMON.NAMES'
2034       include 'COMMON.TIME1'
2035       double precision emgdt(MAXRES6),
2036      & pterm,vterm,rho,rhoc,vsig,
2037      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2038      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2039      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2040       logical lprn /.false./
2041       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2042       double precision ktm
2043 #ifdef MPI
2044       tt0 = MPI_Wtime()
2045 #else
2046       tt0 = tcpu()
2047 #endif
2048 c
2049 c AL 8/17/04 Code adapted from tinker
2050 c
2051 c Get the frictional and random terms for stochastic dynamics in the
2052 c eigenspace of mass-scaled UNRES friction matrix
2053 c
2054       do i = 1, dimen
2055             gdt = fricgam(i) * d_time
2056 c
2057 c Stochastic dynamics reduces to simple MD for zero friction
2058 c
2059             if (gdt .le. zero) then
2060                pfric_vec(i) = 1.0d0
2061                vfric_vec(i) = d_time
2062                afric_vec(i) = 0.5d0 * d_time * d_time
2063                prand_vec(i) = 0.0d0
2064                vrand_vec1(i) = 0.0d0
2065                vrand_vec2(i) = 0.0d0
2066 c
2067 c Analytical expressions when friction coefficient is large
2068 c
2069             else
2070                if (gdt .ge. gdt_radius) then
2071                   egdt = dexp(-gdt)
2072                   pfric_vec(i) = egdt
2073                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2074                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2075                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2076                   vterm = 1.0d0 - egdt**2
2077                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2078 c
2079 c Use series expansions when friction coefficient is small
2080 c
2081                else
2082                   gdt2 = gdt * gdt
2083                   gdt3 = gdt * gdt2
2084                   gdt4 = gdt2 * gdt2
2085                   gdt5 = gdt2 * gdt3
2086                   gdt6 = gdt3 * gdt3
2087                   gdt7 = gdt3 * gdt4
2088                   gdt8 = gdt4 * gdt4
2089                   gdt9 = gdt4 * gdt5
2090                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2091      &                          - gdt5/120.0d0 + gdt6/720.0d0
2092      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2093      &                          - gdt9/362880.0d0) / fricgam(i)**2
2094                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2095                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2096                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2097      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2098      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2099      &                       + 127.0d0*gdt9/90720.0d0
2100                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2101      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2102      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2103      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2104                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2105      &                       - 17.0d0*gdt2/1280.0d0
2106      &                       + 17.0d0*gdt3/6144.0d0
2107      &                       + 40967.0d0*gdt4/34406400.0d0
2108      &                       - 57203.0d0*gdt5/275251200.0d0
2109      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2110                end if
2111 c
2112 c Compute the scaling factors of random terms for the nonzero friction case
2113 c
2114                ktm = 0.5d0*d_time/fricgam(i)
2115                psig = dsqrt(ktm*pterm) / fricgam(i)
2116                vsig = dsqrt(ktm*vterm)
2117                rhoc = dsqrt(1.0d0 - rho*rho)
2118                prand_vec(i) = psig
2119                vrand_vec1(i) = vsig * rho
2120                vrand_vec2(i) = vsig * rhoc
2121             end if
2122       end do
2123       if (lprn) then
2124       write (iout,*)
2125      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2126      &  " vrand_vec2"
2127       do i=1,dimen
2128         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2129      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2130       enddo
2131       endif
2132 c
2133 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2134 c
2135 #ifndef   LANG0
2136       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2137       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2138       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2139       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2140       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2141       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2142 #endif
2143 #ifdef MPI
2144       t_sdsetup=t_sdsetup+MPI_Wtime()
2145 #else
2146       t_sdsetup=t_sdsetup+tcpu()-tt0
2147 #endif
2148       return
2149       end
2150 c-------------------------------------------------------------
2151       subroutine eigtransf1(n,ndim,ab,d,c)
2152       implicit none
2153       integer n,ndim
2154       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2155       integer i,j,k
2156       do i=1,n
2157         do j=1,n
2158           c(i,j)=0.0d0
2159           do k=1,n
2160             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2161           enddo
2162         enddo
2163       enddo
2164       return
2165       end
2166 c-------------------------------------------------------------
2167       subroutine eigtransf(n,ndim,a,b,d,c)
2168       implicit none
2169       integer n,ndim
2170       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2171       integer i,j,k
2172       do i=1,n
2173         do j=1,n
2174           c(i,j)=0.0d0
2175           do k=1,n
2176             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2177           enddo
2178         enddo
2179       enddo
2180       return
2181       end
2182 c-------------------------------------------------------------
2183       subroutine sd_verlet1
2184 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2185       implicit real*8 (a-h,o-z)
2186       include 'DIMENSIONS'
2187       include 'COMMON.CONTROL'
2188       include 'COMMON.VAR'
2189       include 'COMMON.MD'
2190 #ifndef LANG0
2191       include 'COMMON.LANGEVIN'
2192 #else
2193       include 'COMMON.LANGEVIN.lang0'
2194 #endif
2195       include 'COMMON.CHAIN'
2196       include 'COMMON.DERIV'
2197       include 'COMMON.GEO'
2198       include 'COMMON.LOCAL'
2199       include 'COMMON.INTERACT'
2200       include 'COMMON.IOUNITS'
2201       include 'COMMON.NAMES'
2202       double precision stochforcvec(MAXRES6)
2203       common /stochcalc/ stochforcvec
2204       logical lprn /.false./
2205
2206 c      write (iout,*) "dc_old"
2207 c      do i=0,nres
2208 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2209 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2210 c      enddo
2211       do j=1,3
2212         dc_work(j)=dc_old(j,0)
2213         d_t_work(j)=d_t_old(j,0)
2214         d_a_work(j)=d_a_old(j,0)
2215       enddo
2216       ind=3
2217       do i=nnt,nct-1
2218         do j=1,3
2219           dc_work(ind+j)=dc_old(j,i)
2220           d_t_work(ind+j)=d_t_old(j,i)
2221           d_a_work(ind+j)=d_a_old(j,i)
2222         enddo
2223         ind=ind+3
2224       enddo
2225       do i=nnt,nct
2226         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2227           do j=1,3
2228             dc_work(ind+j)=dc_old(j,i+nres)
2229             d_t_work(ind+j)=d_t_old(j,i+nres)
2230             d_a_work(ind+j)=d_a_old(j,i+nres)
2231           enddo
2232           ind=ind+3
2233         endif
2234       enddo
2235 #ifndef LANG0
2236       if (lprn) then
2237       write (iout,*)
2238      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2239      &  " vrand_mat2"
2240       do i=1,dimen
2241         do j=1,dimen
2242           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2243      &      vfric_mat(i,j),afric_mat(i,j),
2244      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2245         enddo
2246       enddo
2247       endif
2248       do i=1,dimen
2249         ddt1=0.0d0
2250         ddt2=0.0d0
2251         do j=1,dimen
2252           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2253      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2254           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2255           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2256         enddo
2257         d_t_work_new(i)=ddt1+0.5d0*ddt2
2258         d_t_work(i)=ddt1+ddt2
2259       enddo
2260 #endif
2261       do j=1,3
2262         dc(j,0)=dc_work(j)
2263         d_t(j,0)=d_t_work(j)
2264       enddo
2265       ind=3
2266       do i=nnt,nct-1
2267         do j=1,3
2268           dc(j,i)=dc_work(ind+j)
2269           d_t(j,i)=d_t_work(ind+j)
2270         enddo
2271         ind=ind+3
2272       enddo
2273       do i=nnt,nct
2274         if (itype(i).ne.10) then
2275           inres=i+nres
2276           do j=1,3
2277             dc(j,inres)=dc_work(ind+j)
2278             d_t(j,inres)=d_t_work(ind+j)
2279           enddo
2280           ind=ind+3
2281         endif
2282       enddo
2283       return
2284       end
2285 c--------------------------------------------------------------------------
2286       subroutine sd_verlet2
2287 c  Calculating the adjusted velocities for accelerations
2288       implicit real*8 (a-h,o-z)
2289       include 'DIMENSIONS'
2290       include 'COMMON.CONTROL'
2291       include 'COMMON.VAR'
2292       include 'COMMON.MD'
2293 #ifndef LANG0
2294       include 'COMMON.LANGEVIN'
2295 #else
2296       include 'COMMON.LANGEVIN.lang0'
2297 #endif
2298       include 'COMMON.CHAIN'
2299       include 'COMMON.DERIV'
2300       include 'COMMON.GEO'
2301       include 'COMMON.LOCAL'
2302       include 'COMMON.INTERACT'
2303       include 'COMMON.IOUNITS'
2304       include 'COMMON.NAMES'
2305       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2306       common /stochcalc/ stochforcvec
2307 c
2308 c Compute the stochastic forces which contribute to velocity change
2309 c
2310       call stochastic_force(stochforcvecV)
2311
2312 #ifndef LANG0
2313       do i=1,dimen
2314         ddt1=0.0d0
2315         ddt2=0.0d0
2316         do j=1,dimen
2317           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2318           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2319      &     vrand_mat2(i,j)*stochforcvecV(j)
2320         enddo
2321         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2322       enddo
2323 #endif
2324       do j=1,3
2325         d_t(j,0)=d_t_work(j)
2326       enddo
2327       ind=3
2328       do i=nnt,nct-1
2329         do j=1,3
2330           d_t(j,i)=d_t_work(ind+j)
2331         enddo
2332         ind=ind+3
2333       enddo
2334       do i=nnt,nct
2335         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2336           inres=i+nres
2337           do j=1,3
2338             d_t(j,inres)=d_t_work(ind+j)
2339           enddo
2340           ind=ind+3
2341         endif
2342       enddo
2343       return
2344       end
2345 c-----------------------------------------------------------
2346       subroutine sd_verlet_ciccotti_setup
2347 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2348 c version
2349       implicit real*8 (a-h,o-z)
2350       include 'DIMENSIONS'
2351 #ifdef MPI
2352       include 'mpif.h'
2353 #endif
2354       include 'COMMON.CONTROL'
2355       include 'COMMON.VAR'
2356       include 'COMMON.MD'
2357 #ifndef LANG0
2358       include 'COMMON.LANGEVIN'
2359 #else
2360       include 'COMMON.LANGEVIN.lang0'
2361 #endif
2362       include 'COMMON.CHAIN'
2363       include 'COMMON.DERIV'
2364       include 'COMMON.GEO'
2365       include 'COMMON.LOCAL'
2366       include 'COMMON.INTERACT'
2367       include 'COMMON.IOUNITS'
2368       include 'COMMON.NAMES'
2369       include 'COMMON.TIME1'
2370       double precision emgdt(MAXRES6),
2371      & pterm,vterm,rho,rhoc,vsig,
2372      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2373      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2374      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2375       logical lprn /.false./
2376       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2377       double precision ktm
2378 #ifdef MPI
2379       tt0 = MPI_Wtime()
2380 #else
2381       tt0 = tcpu()
2382 #endif
2383 c
2384 c AL 8/17/04 Code adapted from tinker
2385 c
2386 c Get the frictional and random terms for stochastic dynamics in the
2387 c eigenspace of mass-scaled UNRES friction matrix
2388 c
2389       do i = 1, dimen
2390             write (iout,*) "i",i," fricgam",fricgam(i)
2391             gdt = fricgam(i) * d_time
2392 c
2393 c Stochastic dynamics reduces to simple MD for zero friction
2394 c
2395             if (gdt .le. zero) then
2396                pfric_vec(i) = 1.0d0
2397                vfric_vec(i) = d_time
2398                afric_vec(i) = 0.5d0*d_time*d_time
2399                prand_vec(i) = afric_vec(i)
2400                vrand_vec2(i) = vfric_vec(i)
2401 c
2402 c Analytical expressions when friction coefficient is large
2403 c
2404             else
2405                egdt = dexp(-gdt)
2406                pfric_vec(i) = egdt
2407                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2408                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2409                prand_vec(i) = afric_vec(i)
2410                vrand_vec2(i) = vfric_vec(i)
2411 c
2412 c Compute the scaling factors of random terms for the nonzero friction case
2413 c
2414 c               ktm = 0.5d0*d_time/fricgam(i)
2415 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2416 c               vsig = dsqrt(ktm*vterm)
2417 c               prand_vec(i) = psig*afric_vec(i)
2418 c               vrand_vec2(i) = vsig*vfric_vec(i)
2419             end if
2420       end do
2421       if (lprn) then
2422       write (iout,*)
2423      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2424      &  " vrand_vec2"
2425       do i=1,dimen
2426         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2427      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2428       enddo
2429       endif
2430 c
2431 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2432 c
2433       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2434       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2435       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2436       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2437       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2438 #ifdef MPI
2439       t_sdsetup=t_sdsetup+MPI_Wtime()
2440 #else
2441       t_sdsetup=t_sdsetup+tcpu()-tt0
2442 #endif
2443       return
2444       end
2445 c-------------------------------------------------------------
2446       subroutine sd_verlet1_ciccotti
2447 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2448       implicit real*8 (a-h,o-z)
2449       include 'DIMENSIONS'
2450 #ifdef MPI
2451       include 'mpif.h'
2452 #endif
2453       include 'COMMON.CONTROL'
2454       include 'COMMON.VAR'
2455       include 'COMMON.MD'
2456 #ifndef LANG0
2457       include 'COMMON.LANGEVIN'
2458 #else
2459       include 'COMMON.LANGEVIN.lang0'
2460 #endif
2461       include 'COMMON.CHAIN'
2462       include 'COMMON.DERIV'
2463       include 'COMMON.GEO'
2464       include 'COMMON.LOCAL'
2465       include 'COMMON.INTERACT'
2466       include 'COMMON.IOUNITS'
2467       include 'COMMON.NAMES'
2468       double precision stochforcvec(MAXRES6)
2469       common /stochcalc/ stochforcvec
2470       logical lprn /.false./
2471
2472 c      write (iout,*) "dc_old"
2473 c      do i=0,nres
2474 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2475 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2476 c      enddo
2477       do j=1,3
2478         dc_work(j)=dc_old(j,0)
2479         d_t_work(j)=d_t_old(j,0)
2480         d_a_work(j)=d_a_old(j,0)
2481       enddo
2482       ind=3
2483       do i=nnt,nct-1
2484         do j=1,3
2485           dc_work(ind+j)=dc_old(j,i)
2486           d_t_work(ind+j)=d_t_old(j,i)
2487           d_a_work(ind+j)=d_a_old(j,i)
2488         enddo
2489         ind=ind+3
2490       enddo
2491       do i=nnt,nct
2492         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2493           do j=1,3
2494             dc_work(ind+j)=dc_old(j,i+nres)
2495             d_t_work(ind+j)=d_t_old(j,i+nres)
2496             d_a_work(ind+j)=d_a_old(j,i+nres)
2497           enddo
2498           ind=ind+3
2499         endif
2500       enddo
2501
2502 #ifndef LANG0
2503       if (lprn) then
2504       write (iout,*)
2505      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2506      &  " vrand_mat2"
2507       do i=1,dimen
2508         do j=1,dimen
2509           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2510      &      vfric_mat(i,j),afric_mat(i,j),
2511      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2512         enddo
2513       enddo
2514       endif
2515       do i=1,dimen
2516         ddt1=0.0d0
2517         ddt2=0.0d0
2518         do j=1,dimen
2519           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2520      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2521           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2522           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2523         enddo
2524         d_t_work_new(i)=ddt1+0.5d0*ddt2
2525         d_t_work(i)=ddt1+ddt2
2526       enddo
2527 #endif
2528       do j=1,3
2529         dc(j,0)=dc_work(j)
2530         d_t(j,0)=d_t_work(j)
2531       enddo
2532       ind=3
2533       do i=nnt,nct-1
2534         do j=1,3
2535           dc(j,i)=dc_work(ind+j)
2536           d_t(j,i)=d_t_work(ind+j)
2537         enddo
2538         ind=ind+3
2539       enddo
2540       do i=nnt,nct
2541         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2542           inres=i+nres
2543           do j=1,3
2544             dc(j,inres)=dc_work(ind+j)
2545             d_t(j,inres)=d_t_work(ind+j)
2546           enddo
2547           ind=ind+3
2548         endif
2549       enddo
2550       return
2551       end
2552 c--------------------------------------------------------------------------
2553       subroutine sd_verlet2_ciccotti
2554 c  Calculating the adjusted velocities for accelerations
2555       implicit real*8 (a-h,o-z)
2556       include 'DIMENSIONS'
2557       include 'COMMON.CONTROL'
2558       include 'COMMON.VAR'
2559       include 'COMMON.MD'
2560 #ifndef LANG0
2561       include 'COMMON.LANGEVIN'
2562 #else
2563       include 'COMMON.LANGEVIN.lang0'
2564 #endif
2565       include 'COMMON.CHAIN'
2566       include 'COMMON.DERIV'
2567       include 'COMMON.GEO'
2568       include 'COMMON.LOCAL'
2569       include 'COMMON.INTERACT'
2570       include 'COMMON.IOUNITS'
2571       include 'COMMON.NAMES'
2572       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2573       common /stochcalc/ stochforcvec
2574 c
2575 c Compute the stochastic forces which contribute to velocity change
2576 c
2577       call stochastic_force(stochforcvecV)
2578 #ifndef LANG0
2579       do i=1,dimen
2580         ddt1=0.0d0
2581         ddt2=0.0d0
2582         do j=1,dimen
2583
2584           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2585 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2586           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2587         enddo
2588         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2589       enddo
2590 #endif
2591       do j=1,3
2592         d_t(j,0)=d_t_work(j)
2593       enddo
2594       ind=3
2595       do i=nnt,nct-1
2596         do j=1,3
2597           d_t(j,i)=d_t_work(ind+j)
2598         enddo
2599         ind=ind+3
2600       enddo
2601       do i=nnt,nct
2602         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2603           inres=i+nres
2604           do j=1,3
2605             d_t(j,inres)=d_t_work(ind+j)
2606           enddo
2607           ind=ind+3
2608         endif
2609       enddo
2610       return
2611       end
2612 #endif