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

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