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