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

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