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

   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       call zerograd
 631       call etotal_short(energia_short)
 632       call cartgrad
 633       call lagrangian
 634       if (ntwe.ne.0) then
 635       if (large.and. mod(itime,ntwe).eq.0) then
 636         write (iout,*) "energia_short",energia_short(0)
 637         write (iout,*) "Accelerations from short-range forces"
 638         do i=0,nres
 639           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 640      &      (d_a(j,i+nres),j=1,3)
 641         enddo
 642       endif
 643       endif
 644 #ifdef MPI
 645         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 646 #else
 647         t_enegrad=t_enegrad+tcpu()-tt0
 648 #endif
 649       do i=0,2*nres
 650         do j=1,3
 651           dc_old(j,i)=dc(j,i)
 652           d_t_old(j,i)=d_t(j,i)
 653           d_a_old(j,i)=d_a(j,i)
 654         enddo
 655       enddo
 656 c 6/30/08 A-MTS: attempt at increasing the split number
 657       do i=0,2*nres
 658         do j=1,3
 659           dc_old0(j,i)=dc_old(j,i)
 660           d_t_old0(j,i)=d_t_old(j,i)
 661           d_a_old0(j,i)=d_a_old(j,i)
 662         enddo
 663       enddo
 664       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 665       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 666 c
 667       scale=.true.
 668       d_time0=d_time
 669       do while (scale)
 670
 671       scale=.false.
 672 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 673 c Split the time step
 674       d_time=d_time0/ntime_split
 675 c Perform the short-range RESPA steps (velocity Verlet increments of
 676 c positions and velocities using short-range forces)
 677 c      write (iout,*) "*********************** RESPA split"
 678       do itsplit=1,ntime_split
 679         if (lang.eq.1) then
 680           call sddir_precalc
 681         else if (lang.eq.2 .or. lang.eq.3) then
 682 #ifndef LANG0
 683           call stochastic_force(stochforcvec)
 684 #else
 685           write (iout,*)
 686      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 687 #ifdef MPI
 688           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 689 #endif
 690           stop
 691 #endif
 692         endif
 693 c First step of the velocity Verlet algorithm
 694         if (lang.eq.2) then
 695 #ifndef LANG0
 696           call sd_verlet1
 697 #endif
 698         else if (lang.eq.3) then
 699 #ifndef LANG0
 700           call sd_verlet1_ciccotti
 701 #endif
 702         else if (lang.eq.1) then
 703           call sddir_verlet1
 704         else
 705           call verlet1
 706         endif
 707 c Build the chain from the newly calculated coordinates
 708         call chainbuild_cart
 709         if (rattle) call rattle1
 710         if (ntwe.ne.0) then
 711         if (large.and. mod(itime,ntwe).eq.0) then
 712           write (iout,*) "***** ITSPLIT",itsplit
 713           write (iout,*) "Cartesian and internal coordinates: step 1"
 714           call pdbout(0.0d0,"cipiszcze",iout)
 715 c          call cartprint
 716           call intout
 717           write (iout,*) "Velocities, step 1"
 718           do i=0,nres
 719             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 720      &        (d_t(j,i+nres),j=1,3)
 721           enddo
 722         endif
 723         endif
 724 #ifdef MPI
 725         tt0 = MPI_Wtime()
 726 #else
 727         tt0 = tcpu()
 728 #endif
 729 c Calculate energy and forces
 730         call zerograd
 731         call etotal_short(energia_short)
 732         call cartgrad
 733 c Get the new accelerations
 734         call lagrangian
 735         if (ntwe.ne.0) then
 736         if (large.and. mod(itime,ntwe).eq.0) then
 737           write (iout,*)"energia_short",energia_short(0)
 738           write (iout,*) "Accelerations from short-range forces"
 739           do i=0,nres
 740             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 741      &        (d_a(j,i+nres),j=1,3)
 742           enddo
 743         endif
 744         endif
 745 c 6/30/08 A-MTS
 746 c Determine maximum acceleration and scale down the timestep if needed
 747         call max_accel
 748         amax=amax/ntime_split**2
 749         call predict_edrift(epdrift)
 750         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 751      &   write (iout,*) "amax",amax," damax",damax,
 752      &   " epdrift",epdrift," epdriftmax",epdriftmax
 753 c Exit loop and try with increased split number if the change of
 754 c acceleration is too big
 755         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 756           if (ntime_split.lt.maxtime_split) then
 757             scale=.true.
 758             ntime_split=ntime_split*2
 759             do i=0,2*nres
 760               do j=1,3
 761                 dc_old(j,i)=dc_old0(j,i)
 762                 d_t_old(j,i)=d_t_old0(j,i)
 763                 d_a_old(j,i)=d_a_old0(j,i)
 764               enddo
 765             enddo
 766             write (iout,*) "acceleration/energy drift too large",amax,
 767      &      epdrift," split increased to ",ntime_split," itime",itime,
 768      &       " itsplit",itsplit
 769             exit
 770           else
 771             write (iout,*)
 772      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 773      &       maxtime_split,
 774      &      " already reached!!! Trying to carry on!"
 775           endif
 776         endif
 777 #ifdef MPI
 778         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 779 #else
 780         t_enegrad=t_enegrad+tcpu()-tt0
 781 #endif
 782 c Second step of the velocity Verlet algorithm
 783         if (lang.eq.2) then
 784 #ifndef LANG0
 785           call sd_verlet2
 786 #endif
 787         else if (lang.eq.3) then
 788 #ifndef LANG0
 789           call sd_verlet2_ciccotti
 790 #endif
 791         else if (lang.eq.1) then
 792           call sddir_verlet2
 793         else
 794           call verlet2
 795         endif
 796         if (rattle) call rattle2
 797 c Backup the coordinates, velocities, and accelerations
 798         do i=0,2*nres
 799           do j=1,3
 800             dc_old(j,i)=dc(j,i)
 801             d_t_old(j,i)=d_t(j,i)
 802             d_a_old(j,i)=d_a(j,i)
 803           enddo
 804         enddo
 805       enddo
 806
 807       enddo ! while scale
 808
 809 c Restore the time step
 810       d_time=d_time0
 811 c Compute long-range forces
 812 #ifdef MPI
 813       tt0 =MPI_Wtime()
 814 #else
 815       tt0 = tcpu()
 816 #endif
 817       call zerograd
 818       call etotal_long(energia_long)
 819       call cartgrad
 820       call lagrangian
 821 #ifdef MPI
 822         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 823 #else
 824         t_enegrad=t_enegrad+tcpu()-tt0
 825 #endif
 826 c Compute accelerations from long-range forces
 827       if (ntwe.ne.0) then
 828       if (large.and. mod(itime,ntwe).eq.0) then
 829         write (iout,*) "energia_long",energia_long(0)
 830         write (iout,*) "Cartesian and internal coordinates: step 2"
 831 c        call cartprint
 832         call pdbout(0.0d0,"cipiszcze",iout)
 833         call intout
 834         write (iout,*) "Accelerations from long-range forces"
 835         do i=0,nres
 836           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 837      &      (d_a(j,i+nres),j=1,3)
 838         enddo
 839         write (iout,*) "Velocities, step 2"
 840         do i=0,nres
 841           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 842      &      (d_t(j,i+nres),j=1,3)
 843         enddo
 844       endif
 845       endif
 846 c Compute the final RESPA step (increment velocities)
 847 c      write (iout,*) "*********************** RESPA fin"
 848       call RESPA_vel
 849 c Compute the complete potential energy
 850       do i=0,n_ene
 851         potEcomp(i)=energia_short(i)+energia_long(i)
 852       enddo
 853       potE=potEcomp(0)-potEcomp(20)
 854 c      potE=energia_short(0)+energia_long(0)
 855       totT=totT+d_time
 856 c Calculate the kinetic and the total energy and the kinetic temperature
 857       call kinetic(EK)
 858       totE=EK+potE
 859 c Couple the system to Berendsen bath if needed
 860       if (tbf .and. lang.eq.0) then
 861         call verlet_bath
 862       endif
 863       kinetic_T=2.0d0/(dimen3*Rb)*EK
 864 c Backup the coordinates, velocities, and accelerations
 865       if (ntwe.ne.0) then
 866       if (mod(itime,ntwe).eq.0 .and. large) then
 867         write (iout,*) "Velocities, end"
 868         do i=0,nres
 869           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 870      &      (d_t(j,i+nres),j=1,3)
 871         enddo
 872       endif
 873       endif
 874       return
 875       end
 876 c---------------------------------------------------------------------
 877       subroutine RESPA_vel
 878 c  First and last RESPA step (incrementing velocities using long-range
 879 c  forces).
 880       implicit real*8 (a-h,o-z)
 881       include 'DIMENSIONS'
 882       include 'COMMON.CONTROL'
 883       include 'COMMON.VAR'
 884       include 'COMMON.MD'
 885       include 'COMMON.CHAIN'
 886       include 'COMMON.DERIV'
 887       include 'COMMON.GEO'
 888       include 'COMMON.LOCAL'
 889       include 'COMMON.INTERACT'
 890       include 'COMMON.IOUNITS'
 891       include 'COMMON.NAMES'
 892       do j=1,3
 893         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
 894       enddo
 895       do i=nnt,nct-1
 896         do j=1,3
 897           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
 898         enddo
 899       enddo
 900       do i=nnt,nct
 901         if (itype(i).ne.10) then
 902           inres=i+nres
 903           do j=1,3
 904             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
 905           enddo
 906         endif
 907       enddo
 908       return
 909       end
 910 c-----------------------------------------------------------------
 911       subroutine verlet1
 912 c Applying velocity Verlet algorithm - step 1 to coordinates
 913       implicit real*8 (a-h,o-z)
 914       include 'DIMENSIONS'
 915       include 'COMMON.CONTROL'
 916       include 'COMMON.VAR'
 917       include 'COMMON.MD'
 918       include 'COMMON.CHAIN'
 919       include 'COMMON.DERIV'
 920       include 'COMMON.GEO'
 921       include 'COMMON.LOCAL'
 922       include 'COMMON.INTERACT'
 923       include 'COMMON.IOUNITS'
 924       include 'COMMON.NAMES'
 925       double precision adt,adt2
 926
 927 #ifdef DEBUG
 928       write (iout,*) "VELVERLET1 START: DC"
 929       do i=0,nres
 930         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
 931      &   (dc(j,i+nres),j=1,3)
 932       enddo
 933 #endif
 934       do j=1,3
 935         adt=d_a_old(j,0)*d_time
 936         adt2=0.5d0*adt
 937         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
 938         d_t_new(j,0)=d_t_old(j,0)+adt2
 939         d_t(j,0)=d_t_old(j,0)+adt
 940       enddo
 941       do i=nnt,nct-1
 942         do j=1,3
 943           adt=d_a_old(j,i)*d_time
 944           adt2=0.5d0*adt
 945           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
 946           d_t_new(j,i)=d_t_old(j,i)+adt2
 947           d_t(j,i)=d_t_old(j,i)+adt
 948         enddo
 949       enddo
 950       do i=nnt,nct
 951         if (itype(i).ne.10) then
 952           inres=i+nres
 953           do j=1,3
 954             adt=d_a_old(j,inres)*d_time
 955             adt2=0.5d0*adt
 956             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
 957             d_t_new(j,inres)=d_t_old(j,inres)+adt2
 958             d_t(j,inres)=d_t_old(j,inres)+adt
 959           enddo
 960         endif
 961       enddo
 962 #ifdef DEBUG
 963       write (iout,*) "VELVERLET1 END: DC"
 964       do i=0,nres
 965         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
 966      &   (dc(j,i+nres),j=1,3)
 967       enddo
 968 #endif
 969       return
 970       end
 971 c---------------------------------------------------------------------
 972       subroutine verlet2
 973 c  Step 2 of the velocity Verlet algorithm: update velocities
 974       implicit real*8 (a-h,o-z)
 975       include 'DIMENSIONS'
 976       include 'COMMON.CONTROL'
 977       include 'COMMON.VAR'
 978       include 'COMMON.MD'
 979       include 'COMMON.CHAIN'
 980       include 'COMMON.DERIV'
 981       include 'COMMON.GEO'
 982       include 'COMMON.LOCAL'
 983       include 'COMMON.INTERACT'
 984       include 'COMMON.IOUNITS'
 985       include 'COMMON.NAMES'
 986       do j=1,3
 987         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
 988       enddo
 989       do i=nnt,nct-1
 990         do j=1,3
 991           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
 992         enddo
 993       enddo
 994       do i=nnt,nct
 995         if (itype(i).ne.10) then
 996           inres=i+nres
 997           do j=1,3
 998             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
 999           enddo
1000         endif
1001       enddo
1002       return
1003       end
1004 c-----------------------------------------------------------------
1005       subroutine sddir_precalc
1006 c Applying velocity Verlet algorithm - step 1 to coordinates
1007       implicit real*8 (a-h,o-z)
1008       include 'DIMENSIONS'
1009       include 'COMMON.CONTROL'
1010       include 'COMMON.VAR'
1011       include 'COMMON.MD'
1012 #ifndef LANG0
1013       include 'COMMON.LANGEVIN'
1014 #else
1015       include 'COMMON.LANGEVIN.lang0'
1016 #endif
1017       include 'COMMON.CHAIN'
1018       include 'COMMON.DERIV'
1019       include 'COMMON.GEO'
1020       include 'COMMON.LOCAL'
1021       include 'COMMON.INTERACT'
1022       include 'COMMON.IOUNITS'
1023       include 'COMMON.NAMES'
1024       double precision stochforcvec(MAXRES6)
1025       common /stochcalc/ stochforcvec
1026 c
1027 c Compute friction and stochastic forces
1028 c
1029       call friction_force
1030       call stochastic_force(stochforcvec)
1031 c
1032 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1033 c forces (d_as_work)
1034 c
1035       call ginv_mult(fric_work, d_af_work)
1036       call ginv_mult(stochforcvec, d_as_work)
1037       return
1038       end
1039 c---------------------------------------------------------------------
1040       subroutine sddir_verlet1
1041 c Applying velocity Verlet algorithm - step 1 to velocities
1042       implicit real*8 (a-h,o-z)
1043       include 'DIMENSIONS'
1044       include 'COMMON.CONTROL'
1045       include 'COMMON.VAR'
1046       include 'COMMON.MD'
1047 #ifndef LANG0
1048       include 'COMMON.LANGEVIN'
1049 #else
1050       include 'COMMON.LANGEVIN.lang0'
1051 #endif
1052       include 'COMMON.CHAIN'
1053       include 'COMMON.DERIV'
1054       include 'COMMON.GEO'
1055       include 'COMMON.LOCAL'
1056       include 'COMMON.INTERACT'
1057       include 'COMMON.IOUNITS'
1058       include 'COMMON.NAMES'
1059 c Revised 3/31/05 AL: correlation between random contributions to
1060 c position and velocity increments included.
1061       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1062       double precision adt,adt2
1063 c
1064 c Add the contribution from BOTH friction and stochastic force to the
1065 c coordinates, but ONLY the contribution from the friction forces to velocities
1066 c
1067       do j=1,3
1068         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1069         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1070         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1071         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1072         d_t(j,0)=d_t_old(j,0)+adt
1073       enddo
1074       ind=3
1075       do i=nnt,nct-1
1076         do j=1,3
1077           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1078           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1079           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1080           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1081           d_t(j,i)=d_t_old(j,i)+adt
1082         enddo
1083         ind=ind+3
1084       enddo
1085       do i=nnt,nct
1086         if (itype(i).ne.10) then
1087           inres=i+nres
1088           do j=1,3
1089             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1090             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1091             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1092             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1093             d_t(j,inres)=d_t_old(j,inres)+adt
1094           enddo
1095           ind=ind+3
1096         endif
1097       enddo
1098       return
1099       end
1100 c---------------------------------------------------------------------
1101       subroutine sddir_verlet2
1102 c  Calculating the adjusted velocities for accelerations
1103       implicit real*8 (a-h,o-z)
1104       include 'DIMENSIONS'
1105       include 'COMMON.CONTROL'
1106       include 'COMMON.VAR'
1107       include 'COMMON.MD'
1108 #ifndef LANG0
1109       include 'COMMON.LANGEVIN'
1110 #else
1111       include 'COMMON.LANGEVIN.lang0'
1112 #endif
1113       include 'COMMON.CHAIN'
1114       include 'COMMON.DERIV'
1115       include 'COMMON.GEO'
1116       include 'COMMON.LOCAL'
1117       include 'COMMON.INTERACT'
1118       include 'COMMON.IOUNITS'
1119       include 'COMMON.NAMES'
1120       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1121       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1122 c Revised 3/31/05 AL: correlation between random contributions to
1123 c position and velocity increments included.
1124 c The correlation coefficients are calculated at low-friction limit.
1125 c Also, friction forces are now not calculated with new velocities.
1126
1127 c      call friction_force
1128       call stochastic_force(stochforcvec)
1129 c
1130 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1131 c forces (d_as_work)
1132 c
1133       call ginv_mult(stochforcvec, d_as_work1)
1134
1135 c
1136 c Update velocities
1137 c
1138       do j=1,3
1139         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1140      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1141       enddo
1142       ind=3
1143       do i=nnt,nct-1
1144         do j=1,3
1145           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1146      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1147         enddo
1148         ind=ind+3
1149       enddo
1150       do i=nnt,nct
1151         if (itype(i).ne.10) then
1152           inres=i+nres
1153           do j=1,3
1154             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1155      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1156      &       +cos60*d_as_work1(ind+j))*d_time
1157           enddo
1158           ind=ind+3
1159         endif
1160       enddo
1161       return
1162       end
1163 c---------------------------------------------------------------------
1164       subroutine max_accel
1165 c
1166 c Find the maximum difference in the accelerations of the the sites
1167 c at the beginning and the end of the time step.
1168 c
1169       implicit real*8 (a-h,o-z)
1170       include 'DIMENSIONS'
1171       include 'COMMON.CONTROL'
1172       include 'COMMON.VAR'
1173       include 'COMMON.MD'
1174       include 'COMMON.CHAIN'
1175       include 'COMMON.DERIV'
1176       include 'COMMON.GEO'
1177       include 'COMMON.LOCAL'
1178       include 'COMMON.INTERACT'
1179       include 'COMMON.IOUNITS'
1180       double precision aux(3),accel(3),accel_old(3),dacc
1181       do j=1,3
1182 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1183          accel_old(j)=d_a_old(j,0)
1184          accel(j)=d_a(j,0)
1185       enddo
1186       amax=0.0d0
1187       do i=nnt,nct
1188 c Backbone
1189         if (i.lt.nct) then
1190 c 7/3/08 changed to asymmetric difference
1191           do j=1,3
1192 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1193             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1194             accel(j)=accel(j)+0.5d0*d_a(j,i)
1195 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1196             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1197               dacc=dabs(accel(j)-accel_old(j))
1198               if (dacc.gt.amax) amax=dacc
1199             endif
1200           enddo
1201         endif
1202       enddo
1203 c Side chains
1204       do j=1,3
1205 c        accel(j)=aux(j)
1206         accel_old(j)=d_a_old(j,0)
1207         accel(j)=d_a(j,0)
1208       enddo
1209       if (nnt.eq.2) then
1210         do j=1,3
1211           accel_old(j)=accel_old(j)+d_a_old(j,1)
1212           accel(j)=accel(j)+d_a(j,1)
1213         enddo
1214       endif
1215       do i=nnt,nct
1216         if (itype(i).ne.10) then
1217           do j=1,3
1218 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1219             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1220             accel(j)=accel(j)+d_a(j,i+nres)
1221           enddo
1222         endif
1223         do j=1,3
1224 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1225           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1226             dacc=dabs(accel(j)-accel_old(j))
1227             if (dacc.gt.amax) amax=dacc
1228           endif
1229         enddo
1230         do j=1,3
1231           accel_old(j)=accel_old(j)+d_a_old(j,i)
1232           accel(j)=accel(j)+d_a(j,i)
1233 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1234         enddo
1235       enddo
1236       return
1237       end
1238 c---------------------------------------------------------------------
1239       subroutine predict_edrift(epdrift)
1240 c
1241 c Predict the drift of the potential energy
1242 c
1243       implicit real*8 (a-h,o-z)
1244       include 'DIMENSIONS'
1245       include 'COMMON.CONTROL'
1246       include 'COMMON.VAR'
1247       include 'COMMON.MD'
1248       include 'COMMON.CHAIN'
1249       include 'COMMON.DERIV'
1250       include 'COMMON.GEO'
1251       include 'COMMON.LOCAL'
1252       include 'COMMON.INTERACT'
1253       include 'COMMON.IOUNITS'
1254       include 'COMMON.MUCA'
1255       double precision epdrift,epdriftij
1256 c Drift of the potential energy
1257       epdrift=0.0d0
1258       do i=nnt,nct
1259 c Backbone
1260         if (i.lt.nct) then
1261           do j=1,3
1262             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1263             if (lmuca) epdriftij=epdriftij*factor
1264 c            write (iout,*) "back",i,j,epdriftij
1265             if (epdriftij.gt.epdrift) epdrift=epdriftij
1266           enddo
1267         endif
1268 c Side chains
1269         if (itype(i).ne.10) then
1270           do j=1,3
1271             epdriftij=
1272      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1273             if (lmuca) epdriftij=epdriftij*factor
1274 c            write (iout,*) "side",i,j,epdriftij
1275             if (epdriftij.gt.epdrift) epdrift=epdriftij
1276           enddo
1277         endif
1278       enddo
1279       epdrift=0.5d0*epdrift*d_time*d_time
1280 c      write (iout,*) "epdrift",epdrift
1281       return
1282       end
1283 c-----------------------------------------------------------------------
1284       subroutine verlet_bath
1285 c
1286 c  Coupling to the thermostat by using the Berendsen algorithm
1287 c
1288       implicit real*8 (a-h,o-z)
1289       include 'DIMENSIONS'
1290       include 'COMMON.CONTROL'
1291       include 'COMMON.VAR'
1292       include 'COMMON.MD'
1293       include 'COMMON.CHAIN'
1294       include 'COMMON.DERIV'
1295       include 'COMMON.GEO'
1296       include 'COMMON.LOCAL'
1297       include 'COMMON.INTERACT'
1298       include 'COMMON.IOUNITS'
1299       include 'COMMON.NAMES'
1300       double precision T_half,fact
1301 c
1302       T_half=2.0d0/(dimen3*Rb)*EK
1303       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1304 c      write(iout,*) "T_half", T_half
1305 c      write(iout,*) "EK", EK
1306 c      write(iout,*) "fact", fact
1307       do j=1,3
1308         d_t(j,0)=fact*d_t(j,0)
1309       enddo
1310       do i=nnt,nct-1
1311         do j=1,3
1312           d_t(j,i)=fact*d_t(j,i)
1313         enddo
1314       enddo
1315       do i=nnt,nct
1316         if (itype(i).ne.10) then
1317           inres=i+nres
1318           do j=1,3
1319             d_t(j,inres)=fact*d_t(j,inres)
1320           enddo
1321         endif
1322       enddo
1323       return
1324       end
1325 c---------------------------------------------------------
1326       subroutine init_MD
1327 c  Set up the initial conditions of a MD simulation
1328       implicit real*8 (a-h,o-z)
1329       include 'DIMENSIONS'
1330 #ifdef MP
1331       include 'mpif.h'
1332       character*16 form
1333       integer IERROR,ERRCODE
1334 #endif
1335       include 'COMMON.SETUP'
1336       include 'COMMON.CONTROL'
1337       include 'COMMON.VAR'
1338       include 'COMMON.MD'
1339 #ifndef LANG0
1340       include 'COMMON.LANGEVIN'
1341 #else
1342       include 'COMMON.LANGEVIN.lang0'
1343 #endif
1344       include 'COMMON.CHAIN'
1345       include 'COMMON.DERIV'
1346       include 'COMMON.GEO'
1347       include 'COMMON.LOCAL'
1348       include 'COMMON.INTERACT'
1349       include 'COMMON.IOUNITS'
1350       include 'COMMON.NAMES'
1351       include 'COMMON.REMD'
1352       real*8 energia_long(0:n_ene),
1353      &  energia_short(0:n_ene),vcm(3),incr(3)
1354       double precision cm(3),L(3),xv,sigv,lowb,highb
1355       double precision varia(maxvar)
1356       character*256 qstr
1357       integer ilen
1358       external ilen
1359       character*50 tytul
1360       logical file_exist
1361       common /gucio/ cm
1362       d_time0=d_time
1363 c      write(iout,*) "d_time", d_time
1364 c Compute the standard deviations of stochastic forces for Langevin dynamics
1365 c if the friction coefficients do not depend on surface area
1366       if (lang.gt.0 .and. .not.surfarea) then
1367         do i=nnt,nct-1
1368           stdforcp(i)=stdfp*dsqrt(gamp)
1369         enddo
1370         do i=nnt,nct
1371           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1372         enddo
1373       endif
1374 c Open the pdb file for snapshotshots
1375 #ifdef MPI
1376       if(mdpdb) then
1377         if (ilen(tmpdir).gt.0)
1378      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1379      &      liczba(:ilen(liczba))//".pdb")
1380         open(ipdb,
1381      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1382      &  //".pdb")
1383       else
1384 #ifdef NOXDR
1385         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1386      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1387      &      liczba(:ilen(liczba))//".x")
1388         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1389      &  //".x"
1390 #else
1391         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1392      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1393      &      liczba(:ilen(liczba))//".cx")
1394         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1395      &  //".cx"
1396 #endif
1397       endif
1398 #else
1399       if(mdpdb) then
1400          if (ilen(tmpdir).gt.0)
1401      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1402          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1403       else
1404          if (ilen(tmpdir).gt.0)
1405      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1406          cartname=prefix(:ilen(prefix))//"_MD.cx"
1407       endif
1408 #endif
1409       if (usampl) then
1410         write (qstr,'(256(1h ))')
1411         ipos=1
1412         do i=1,nfrag
1413           iq = qinfrag(i,iset)*10
1414           iw = wfrag(i,iset)/100
1415           if (iw.gt.0) then
1416             if(me.eq.king.or..not.out1file)
1417      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1418             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1419             ipos=ipos+7
1420           endif
1421         enddo
1422         do i=1,npair
1423           iq = qinpair(i,iset)*10
1424           iw = wpair(i,iset)/100
1425           if (iw.gt.0) then
1426             if(me.eq.king.or..not.out1file)
1427      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1428             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1429             ipos=ipos+7
1430           endif
1431         enddo
1432 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1433 #ifdef NOXDR
1434 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1435 #else
1436 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1437 #endif
1438 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1439       endif
1440       icg=1
1441       if (rest) then
1442        if (restart1file) then
1443          if (me.eq.king)
1444      &     inquire(file=mremd_rst_name,exist=file_exist)
1445            write (*,*) me," Before broadcast: file_exist",file_exist
1446          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1447      &          IERR)
1448          write (*,*) me," After broadcast: file_exist",file_exist
1449 c        inquire(file=mremd_rst_name,exist=file_exist)
1450         if(me.eq.king.or..not.out1file)
1451      &   write(iout,*) "Initial state read by master and distributed"
1452        else
1453          if (ilen(tmpdir).gt.0)
1454      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1455      &      //liczba(:ilen(liczba))//'.rst')
1456         inquire(file=rest2name,exist=file_exist)
1457        endif
1458        if(file_exist) then
1459          if(.not.restart1file) then
1460            if(me.eq.king.or..not.out1file)
1461      &      write(iout,*) "Initial state will be read from file ",
1462      &      rest2name(:ilen(rest2name))
1463            call readrst
1464          endif
1465          call rescale_weights(t_bath)
1466        else
1467         if(me.eq.king.or..not.out1file)then
1468          if (restart1file) then
1469           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1470      &       " does not exist"
1471          else
1472           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1473      &       " does not exist"
1474          endif
1475          write(iout,*) "Initial velocities randomly generated"
1476         endif
1477         call random_vel
1478         totT=0.0d0
1479        endif
1480       else
1481 c Generate initial velocities
1482         if(me.eq.king.or..not.out1file)
1483      &   write(iout,*) "Initial velocities randomly generated"
1484         call random_vel
1485         totT=0.0d0
1486       endif
1487 c      rest2name = prefix(:ilen(prefix))//'.rst'
1488       if(me.eq.king.or..not.out1file)then
1489        write(iout,*) "Initial backbone velocities"
1490        do i=nnt,nct-1
1491         write(iout,*) (d_t(j,i),j=1,3)
1492        enddo
1493        write(iout,*) "Initial side-chain velocities"
1494        do i=nnt,nct
1495         write(iout,*) (d_t(j,i+nres),j=1,3)
1496        enddo
1497 c  Zeroing the total angular momentum of the system
1498        write(iout,*) "Calling the zero-angular
1499      & momentum subroutine"
1500       endif
1501       call inertia_tensor
1502 c  Getting the potential energy and forces and velocities and accelerations
1503       call vcm_vel(vcm)
1504 c      write (iout,*) "velocity of the center of the mass:"
1505 c      write (iout,*) (vcm(j),j=1,3)
1506       do j=1,3
1507         d_t(j,0)=d_t(j,0)-vcm(j)
1508       enddo
1509 c Removing the velocity of the center of mass
1510       call vcm_vel(vcm)
1511       if(me.eq.king.or..not.out1file)then
1512        write (iout,*) "vcm right after adjustment:"
1513        write (iout,*) (vcm(j),j=1,3)
1514       endif
1515       if (.not.rest) then
1516          call chainbuild
1517          if(iranconf.ne.0) then
1518           if (overlapsc) then
1519            print *, 'Calling OVERLAP_SC'
1520            call overlap_sc(fail)
1521           endif
1522
1523           if (searchsc) then
1524            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1525            print *,'SC_move',nft_sc,etot
1526            if(me.eq.king.or..not.out1file)
1527      &      write(iout,*) 'SC_move',nft_sc,etot
1528           endif
1529
1530           if(dccart)then
1531            print *, 'Calling MINIM_DC'
1532            call minim_dc(etot,iretcode,nfun)
1533           else
1534            call geom_to_var(nvar,varia)
1535            print *,'Calling MINIMIZE.'
1536            call minimize(etot,varia,iretcode,nfun)
1537            call var_to_geom(nvar,varia)
1538           endif
1539           if(me.eq.king.or..not.out1file)
1540      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1541          endif
1542       endif
1543       call chainbuild_cart
1544       call kinetic(EK)
1545       if (tbf) then
1546         call verlet_bath(EK)
1547       endif
1548       kinetic_T=2.0d0/(dimen3*Rb)*EK
1549       if(me.eq.king.or..not.out1file)then
1550        call cartprint
1551        call intout
1552       endif
1553       call zerograd
1554       call etotal(potEcomp)
1555       potE=potEcomp(0)
1556       call cartgrad
1557       call lagrangian
1558       call max_accel
1559       if (amax*d_time .gt. dvmax) then
1560         d_time=d_time*dvmax/amax
1561         if(me.eq.king.or..not.out1file) write (iout,*)
1562      &   "Time step reduced to",d_time,
1563      &   " because of too large initial acceleration."
1564       endif
1565       if(me.eq.king.or..not.out1file)then
1566        write(iout,*) "Potential energy and its components"
1567        write(iout,*) (potEcomp(i),i=0,n_ene)
1568       endif
1569       potE=potEcomp(0)-potEcomp(20)
1570       totE=EK+potE
1571       itime=0
1572       if (ntwe.ne.0) call statout(itime)
1573       if(me.eq.king.or..not.out1file)
1574      &  write (iout,*) "Initial:",
1575      &   " Kinetic energy",EK," potential energy",potE,
1576      &   " total energy",totE," maximum acceleration ",
1577      &   amax
1578       if (large) then
1579         write (iout,*) "Initial coordinates"
1580         do i=1,nres
1581           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1582      &    (c(j,i+nres),j=1,3)
1583         enddo
1584         write (iout,*) "Initial dC"
1585         do i=0,nres
1586           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1587      &    (dc(j,i+nres),j=1,3)
1588         enddo
1589         write (iout,*) "Initial velocities"
1590         do i=0,nres
1591           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1592      &    (d_t(j,i+nres),j=1,3)
1593         enddo
1594         write (iout,*) "Initial accelerations"
1595         do i=0,nres
1596 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1597           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1598      &    (d_a(j,i+nres),j=1,3)
1599         enddo
1600       endif
1601       do i=0,2*nres
1602         do j=1,3
1603           dc_old(j,i)=dc(j,i)
1604           d_t_old(j,i)=d_t(j,i)
1605           d_a_old(j,i)=d_a(j,i)
1606         enddo
1607 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1608       enddo
1609       if (RESPA) then
1610 #ifdef MPI
1611       tt0 =MPI_Wtime()
1612 #else
1613       tt0 = tcpu()
1614 #endif
1615         call zerograd
1616         call etotal_long(energia_long)
1617         if(.not.out1file .and. large) then
1618           write (iout,*) "energia_long",energia_long(0)
1619           write (iout,*) "Initial slow-force accelerations"
1620           do i=0,nres
1621             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1622      &      (d_a(j,i+nres),j=1,3)
1623           enddo
1624         endif
1625         call cartgrad
1626         call lagrangian
1627 #ifdef MPI
1628         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1629 #else
1630         t_enegrad=t_enegrad+tcpu()-tt0
1631 #endif
1632 c        call etotal_short(energia_short)
1633 c        write (iout,*) "energia_long",energia_long(0),
1634 c     &   " energia_short",energia_short(0),
1635 c     &   " total",energia_long(0)+energia_short(0)
1636       endif
1637       return
1638       end
1639 c-----------------------------------------------------------
1640       subroutine random_vel
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643       include 'COMMON.CONTROL'
1644       include 'COMMON.VAR'
1645       include 'COMMON.MD'
1646 #ifndef LANG0
1647       include 'COMMON.LANGEVIN'
1648 #else
1649       include 'COMMON.LANGEVIN.lang0'
1650 #endif
1651       include 'COMMON.CHAIN'
1652       include 'COMMON.DERIV'
1653       include 'COMMON.GEO'
1654       include 'COMMON.LOCAL'
1655       include 'COMMON.INTERACT'
1656       include 'COMMON.IOUNITS'
1657       include 'COMMON.NAMES'
1658       include 'COMMON.TIME1'
1659       double precision xv,sigv,lowb,highb
1660 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1661 c First generate velocities in the eigenspace of the G matrix
1662 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1663 c      call flush(iout)
1664       xv=0.0d0
1665       ii=0
1666       do i=1,dimen
1667         do k=1,3
1668           ii=ii+1
1669           sigv=dsqrt((Rb*t_bath)/geigen(i))
1670           lowb=-5*sigv
1671           highb=5*sigv
1672           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1673 c          write (iout,*) "ii",ii," d_t_work_new",d_t_work_new(ii)
1674         enddo
1675       enddo
1676 c diagnostics
1677 c      Ek1=0.0d0
1678 c      ii=0
1679 c      do i=1,dimen
1680 c        do k=1,3
1681 c          ii=ii+1
1682 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1683 c        enddo
1684 c      enddo
1685 c      write (iout,*) "Ek from eigenvectors",Ek1
1686 c end diagnostics
1687 c Transform velocities to UNRES coordinate space
1688       do k=0,2
1689         do i=1,dimen
1690           ind=(i-1)*3+k+1
1691           d_t_work(ind)=0.0d0
1692           do j=1,dimen
1693             d_t_work(ind)=d_t_work(ind)
1694      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1695           enddo
1696 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1697 c          call flush(iout)
1698         enddo
1699       enddo
1700 c Transfer to the d_t vector
1701       do j=1,3
1702         d_t(j,0)=d_t_work(j)
1703       enddo
1704       ind=3
1705       do i=nnt,nct-1
1706         do j=1,3
1707           ind=ind+1
1708           d_t(j,i)=d_t_work(ind)
1709         enddo
1710       enddo
1711       do i=nnt,nct
1712         if (itype(i).ne.10) then
1713           do j=1,3
1714             ind=ind+1
1715             d_t(j,i+nres)=d_t_work(ind)
1716           enddo
1717         endif
1718       enddo
1719 c      call kinetic(EK)
1720 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1721 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1722 c      call flush(iout)
1723       return
1724       end
1725 #ifndef LANG0
1726 c-----------------------------------------------------------
1727       subroutine sd_verlet_p_setup
1728 c Sets up the parameters of stochastic Verlet algorithm
1729       implicit real*8 (a-h,o-z)
1730       include 'DIMENSIONS'
1731 #ifdef MPI
1732       include 'mpif.h'
1733 #endif
1734       include 'COMMON.CONTROL'
1735       include 'COMMON.VAR'
1736       include 'COMMON.MD'
1737 #ifndef LANG0
1738       include 'COMMON.LANGEVIN'
1739 #else
1740       include 'COMMON.LANGEVIN.lang0'
1741 #endif
1742       include 'COMMON.CHAIN'
1743       include 'COMMON.DERIV'
1744       include 'COMMON.GEO'
1745       include 'COMMON.LOCAL'
1746       include 'COMMON.INTERACT'
1747       include 'COMMON.IOUNITS'
1748       include 'COMMON.NAMES'
1749       include 'COMMON.TIME1'
1750       double precision emgdt(MAXRES6),
1751      & pterm,vterm,rho,rhoc,vsig,
1752      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1753      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1754      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1755       logical lprn /.false./
1756       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1757       double precision ktm
1758 #ifdef MPI
1759       tt0 = MPI_Wtime()
1760 #else
1761       tt0 = tcpu()
1762 #endif
1763 c
1764 c AL 8/17/04 Code adapted from tinker
1765 c
1766 c Get the frictional and random terms for stochastic dynamics in the
1767 c eigenspace of mass-scaled UNRES friction matrix
1768 c
1769       do i = 1, dimen
1770             gdt = fricgam(i) * d_time
1771 c
1772 c Stochastic dynamics reduces to simple MD for zero friction
1773 c
1774             if (gdt .le. zero) then
1775                pfric_vec(i) = 1.0d0
1776                vfric_vec(i) = d_time
1777                afric_vec(i) = 0.5d0 * d_time * d_time
1778                prand_vec(i) = 0.0d0
1779                vrand_vec1(i) = 0.0d0
1780                vrand_vec2(i) = 0.0d0
1781 c
1782 c Analytical expressions when friction coefficient is large
1783 c
1784             else
1785                if (gdt .ge. gdt_radius) then
1786                   egdt = dexp(-gdt)
1787                   pfric_vec(i) = egdt
1788                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1789                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1790                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1791                   vterm = 1.0d0 - egdt**2
1792                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1793 c
1794 c Use series expansions when friction coefficient is small
1795 c
1796                else
1797                   gdt2 = gdt * gdt
1798                   gdt3 = gdt * gdt2
1799                   gdt4 = gdt2 * gdt2
1800                   gdt5 = gdt2 * gdt3
1801                   gdt6 = gdt3 * gdt3
1802                   gdt7 = gdt3 * gdt4
1803                   gdt8 = gdt4 * gdt4
1804                   gdt9 = gdt4 * gdt5
1805                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1806      &                          - gdt5/120.0d0 + gdt6/720.0d0
1807      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
1808      &                          - gdt9/362880.0d0) / fricgam(i)**2
1809                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1810                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1811                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1812      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1813      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1814      &                       + 127.0d0*gdt9/90720.0d0
1815                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1816      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1817      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1818      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1819                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1820      &                       - 17.0d0*gdt2/1280.0d0
1821      &                       + 17.0d0*gdt3/6144.0d0
1822      &                       + 40967.0d0*gdt4/34406400.0d0
1823      &                       - 57203.0d0*gdt5/275251200.0d0
1824      &                       - 1429487.0d0*gdt6/13212057600.0d0)
1825                end if
1826 c
1827 c Compute the scaling factors of random terms for the nonzero friction case
1828 c
1829                ktm = 0.5d0*d_time/fricgam(i)
1830                psig = dsqrt(ktm*pterm) / fricgam(i)
1831                vsig = dsqrt(ktm*vterm)
1832                rhoc = dsqrt(1.0d0 - rho*rho)
1833                prand_vec(i) = psig
1834                vrand_vec1(i) = vsig * rho
1835                vrand_vec2(i) = vsig * rhoc
1836             end if
1837       end do
1838       if (lprn) then
1839       write (iout,*)
1840      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1841      &  " vrand_vec2"
1842       do i=1,dimen
1843         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1844      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1845       enddo
1846       endif
1847 c
1848 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1849 c
1850 #ifndef   LANG0
1851       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1852       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1853       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1854       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1855       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1856       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1857 #endif
1858 #ifdef MPI
1859       t_sdsetup=t_sdsetup+MPI_Wtime()
1860 #else
1861       t_sdsetup=t_sdsetup+tcpu()-tt0
1862 #endif
1863       return
1864       end
1865 c-------------------------------------------------------------
1866       subroutine eigtransf1(n,ndim,ab,d,c)
1867       implicit none
1868       integer n,ndim
1869       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1870       integer i,j,k
1871       do i=1,n
1872         do j=1,n
1873           c(i,j)=0.0d0
1874           do k=1,n
1875             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1876           enddo
1877         enddo
1878       enddo
1879       return
1880       end
1881 c-------------------------------------------------------------
1882       subroutine eigtransf(n,ndim,a,b,d,c)
1883       implicit none
1884       integer n,ndim
1885       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
1886       integer i,j,k
1887       do i=1,n
1888         do j=1,n
1889           c(i,j)=0.0d0
1890           do k=1,n
1891             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
1892           enddo
1893         enddo
1894       enddo
1895       return
1896       end
1897 c-------------------------------------------------------------
1898       subroutine sd_verlet1
1899 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
1900       implicit real*8 (a-h,o-z)
1901       include 'DIMENSIONS'
1902       include 'COMMON.CONTROL'
1903       include 'COMMON.VAR'
1904       include 'COMMON.MD'
1905 #ifndef LANG0
1906       include 'COMMON.LANGEVIN'
1907 #else
1908       include 'COMMON.LANGEVIN.lang0'
1909 #endif
1910       include 'COMMON.CHAIN'
1911       include 'COMMON.DERIV'
1912       include 'COMMON.GEO'
1913       include 'COMMON.LOCAL'
1914       include 'COMMON.INTERACT'
1915       include 'COMMON.IOUNITS'
1916       include 'COMMON.NAMES'
1917       double precision stochforcvec(MAXRES6)
1918       common /stochcalc/ stochforcvec
1919       logical lprn /.false./
1920
1921 c      write (iout,*) "dc_old"
1922 c      do i=0,nres
1923 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
1924 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
1925 c      enddo
1926       do j=1,3
1927         dc_work(j)=dc_old(j,0)
1928         d_t_work(j)=d_t_old(j,0)
1929         d_a_work(j)=d_a_old(j,0)
1930       enddo
1931       ind=3
1932       do i=nnt,nct-1
1933         do j=1,3
1934           dc_work(ind+j)=dc_old(j,i)
1935           d_t_work(ind+j)=d_t_old(j,i)
1936           d_a_work(ind+j)=d_a_old(j,i)
1937         enddo
1938         ind=ind+3
1939       enddo
1940       do i=nnt,nct
1941         if (itype(i).ne.10) then
1942           do j=1,3
1943             dc_work(ind+j)=dc_old(j,i+nres)
1944             d_t_work(ind+j)=d_t_old(j,i+nres)
1945             d_a_work(ind+j)=d_a_old(j,i+nres)
1946           enddo
1947           ind=ind+3
1948         endif
1949       enddo
1950 #ifndef LANG0
1951       if (lprn) then
1952       write (iout,*)
1953      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
1954      &  " vrand_mat2"
1955       do i=1,dimen
1956         do j=1,dimen
1957           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
1958      &      vfric_mat(i,j),afric_mat(i,j),
1959      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
1960         enddo
1961       enddo
1962       endif
1963       do i=1,dimen
1964         ddt1=0.0d0
1965         ddt2=0.0d0
1966         do j=1,dimen
1967           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
1968      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
1969           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
1970           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
1971         enddo
1972         d_t_work_new(i)=ddt1+0.5d0*ddt2
1973         d_t_work(i)=ddt1+ddt2
1974       enddo
1975 #endif
1976       do j=1,3
1977         dc(j,0)=dc_work(j)
1978         d_t(j,0)=d_t_work(j)
1979       enddo
1980       ind=3
1981       do i=nnt,nct-1
1982         do j=1,3
1983           dc(j,i)=dc_work(ind+j)
1984           d_t(j,i)=d_t_work(ind+j)
1985         enddo
1986         ind=ind+3
1987       enddo
1988       do i=nnt,nct
1989         if (itype(i).ne.10) then
1990           inres=i+nres
1991           do j=1,3
1992             dc(j,inres)=dc_work(ind+j)
1993             d_t(j,inres)=d_t_work(ind+j)
1994           enddo
1995           ind=ind+3
1996         endif
1997       enddo
1998       return
1999       end
2000 c--------------------------------------------------------------------------
2001       subroutine sd_verlet2
2002 c  Calculating the adjusted velocities for accelerations
2003       implicit real*8 (a-h,o-z)
2004       include 'DIMENSIONS'
2005       include 'COMMON.CONTROL'
2006       include 'COMMON.VAR'
2007       include 'COMMON.MD'
2008 #ifndef LANG0
2009       include 'COMMON.LANGEVIN'
2010 #else
2011       include 'COMMON.LANGEVIN.lang0'
2012 #endif
2013       include 'COMMON.CHAIN'
2014       include 'COMMON.DERIV'
2015       include 'COMMON.GEO'
2016       include 'COMMON.LOCAL'
2017       include 'COMMON.INTERACT'
2018       include 'COMMON.IOUNITS'
2019       include 'COMMON.NAMES'
2020       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2021       common /stochcalc/ stochforcvec
2022 c
2023 c Compute the stochastic forces which contribute to velocity change
2024 c
2025       call stochastic_force(stochforcvecV)
2026
2027 #ifndef LANG0
2028       do i=1,dimen
2029         ddt1=0.0d0
2030         ddt2=0.0d0
2031         do j=1,dimen
2032           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2033           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2034      &     vrand_mat2(i,j)*stochforcvecV(j)
2035         enddo
2036         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2037       enddo
2038 #endif
2039       do j=1,3
2040         d_t(j,0)=d_t_work(j)
2041       enddo
2042       ind=3
2043       do i=nnt,nct-1
2044         do j=1,3
2045           d_t(j,i)=d_t_work(ind+j)
2046         enddo
2047         ind=ind+3
2048       enddo
2049       do i=nnt,nct
2050         if (itype(i).ne.10) then
2051           inres=i+nres
2052           do j=1,3
2053             d_t(j,inres)=d_t_work(ind+j)
2054           enddo
2055           ind=ind+3
2056         endif
2057       enddo
2058       return
2059       end
2060 c-----------------------------------------------------------
2061       subroutine sd_verlet_ciccotti_setup
2062 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2063 c version
2064       implicit real*8 (a-h,o-z)
2065       include 'DIMENSIONS'
2066 #ifdef MPI
2067       include 'mpif.h'
2068 #endif
2069       include 'COMMON.CONTROL'
2070       include 'COMMON.VAR'
2071       include 'COMMON.MD'
2072 #ifndef LANG0
2073       include 'COMMON.LANGEVIN'
2074 #else
2075       include 'COMMON.LANGEVIN.lang0'
2076 #endif
2077       include 'COMMON.CHAIN'
2078       include 'COMMON.DERIV'
2079       include 'COMMON.GEO'
2080       include 'COMMON.LOCAL'
2081       include 'COMMON.INTERACT'
2082       include 'COMMON.IOUNITS'
2083       include 'COMMON.NAMES'
2084       include 'COMMON.TIME1'
2085       double precision emgdt(MAXRES6),
2086      & pterm,vterm,rho,rhoc,vsig,
2087      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2088      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2089      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2090       logical lprn /.false./
2091       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2092       double precision ktm
2093 #ifdef MPI
2094       tt0 = MPI_Wtime()
2095 #else
2096       tt0 = tcpu()
2097 #endif
2098 c
2099 c AL 8/17/04 Code adapted from tinker
2100 c
2101 c Get the frictional and random terms for stochastic dynamics in the
2102 c eigenspace of mass-scaled UNRES friction matrix
2103 c
2104       do i = 1, dimen
2105             write (iout,*) "i",i," fricgam",fricgam(i)
2106             gdt = fricgam(i) * d_time
2107 c
2108 c Stochastic dynamics reduces to simple MD for zero friction
2109 c
2110             if (gdt .le. zero) then
2111                pfric_vec(i) = 1.0d0
2112                vfric_vec(i) = d_time
2113                afric_vec(i) = 0.5d0*d_time*d_time
2114                prand_vec(i) = afric_vec(i)
2115                vrand_vec2(i) = vfric_vec(i)
2116 c
2117 c Analytical expressions when friction coefficient is large
2118 c
2119             else
2120                egdt = dexp(-gdt)
2121                pfric_vec(i) = egdt
2122                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2123                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2124                prand_vec(i) = afric_vec(i)
2125                vrand_vec2(i) = vfric_vec(i)
2126 c
2127 c Compute the scaling factors of random terms for the nonzero friction case
2128 c
2129 c               ktm = 0.5d0*d_time/fricgam(i)
2130 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2131 c               vsig = dsqrt(ktm*vterm)
2132 c               prand_vec(i) = psig*afric_vec(i)
2133 c               vrand_vec2(i) = vsig*vfric_vec(i)
2134             end if
2135       end do
2136       if (lprn) then
2137       write (iout,*)
2138      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2139      &  " vrand_vec2"
2140       do i=1,dimen
2141         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2142      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2143       enddo
2144       endif
2145 c
2146 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2147 c
2148       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2149       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2150       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2151       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2152       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2153 #ifdef MPI
2154       t_sdsetup=t_sdsetup+MPI_Wtime()
2155 #else
2156       t_sdsetup=t_sdsetup+tcpu()-tt0
2157 #endif
2158       return
2159       end
2160 c-------------------------------------------------------------
2161       subroutine sd_verlet1_ciccotti
2162 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2163       implicit real*8 (a-h,o-z)
2164       include 'DIMENSIONS'
2165 #ifdef MPI
2166       include 'mpif.h'
2167 #endif
2168       include 'COMMON.CONTROL'
2169       include 'COMMON.VAR'
2170       include 'COMMON.MD'
2171 #ifndef LANG0
2172       include 'COMMON.LANGEVIN'
2173 #else
2174       include 'COMMON.LANGEVIN.lang0'
2175 #endif
2176       include 'COMMON.CHAIN'
2177       include 'COMMON.DERIV'
2178       include 'COMMON.GEO'
2179       include 'COMMON.LOCAL'
2180       include 'COMMON.INTERACT'
2181       include 'COMMON.IOUNITS'
2182       include 'COMMON.NAMES'
2183       double precision stochforcvec(MAXRES6)
2184       common /stochcalc/ stochforcvec
2185       logical lprn /.false./
2186
2187 c      write (iout,*) "dc_old"
2188 c      do i=0,nres
2189 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2190 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2191 c      enddo
2192       do j=1,3
2193         dc_work(j)=dc_old(j,0)
2194         d_t_work(j)=d_t_old(j,0)
2195         d_a_work(j)=d_a_old(j,0)
2196       enddo
2197       ind=3
2198       do i=nnt,nct-1
2199         do j=1,3
2200           dc_work(ind+j)=dc_old(j,i)
2201           d_t_work(ind+j)=d_t_old(j,i)
2202           d_a_work(ind+j)=d_a_old(j,i)
2203         enddo
2204         ind=ind+3
2205       enddo
2206       do i=nnt,nct
2207         if (itype(i).ne.10) then
2208           do j=1,3
2209             dc_work(ind+j)=dc_old(j,i+nres)
2210             d_t_work(ind+j)=d_t_old(j,i+nres)
2211             d_a_work(ind+j)=d_a_old(j,i+nres)
2212           enddo
2213           ind=ind+3
2214         endif
2215       enddo
2216
2217 #ifndef LANG0
2218       if (lprn) then
2219       write (iout,*)
2220      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2221      &  " vrand_mat2"
2222       do i=1,dimen
2223         do j=1,dimen
2224           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2225      &      vfric_mat(i,j),afric_mat(i,j),
2226      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2227         enddo
2228       enddo
2229       endif
2230       do i=1,dimen
2231         ddt1=0.0d0
2232         ddt2=0.0d0
2233         do j=1,dimen
2234           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2235      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2236           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2237           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2238         enddo
2239         d_t_work_new(i)=ddt1+0.5d0*ddt2
2240         d_t_work(i)=ddt1+ddt2
2241       enddo
2242 #endif
2243       do j=1,3
2244         dc(j,0)=dc_work(j)
2245         d_t(j,0)=d_t_work(j)
2246       enddo
2247       ind=3
2248       do i=nnt,nct-1
2249         do j=1,3
2250           dc(j,i)=dc_work(ind+j)
2251           d_t(j,i)=d_t_work(ind+j)
2252         enddo
2253         ind=ind+3
2254       enddo
2255       do i=nnt,nct
2256         if (itype(i).ne.10) then
2257           inres=i+nres
2258           do j=1,3
2259             dc(j,inres)=dc_work(ind+j)
2260             d_t(j,inres)=d_t_work(ind+j)
2261           enddo
2262           ind=ind+3
2263         endif
2264       enddo
2265       return
2266       end
2267 c--------------------------------------------------------------------------
2268       subroutine sd_verlet2_ciccotti
2269 c  Calculating the adjusted velocities for accelerations
2270       implicit real*8 (a-h,o-z)
2271       include 'DIMENSIONS'
2272       include 'COMMON.CONTROL'
2273       include 'COMMON.VAR'
2274       include 'COMMON.MD'
2275 #ifndef LANG0
2276       include 'COMMON.LANGEVIN'
2277 #else
2278       include 'COMMON.LANGEVIN.lang0'
2279 #endif
2280       include 'COMMON.CHAIN'
2281       include 'COMMON.DERIV'
2282       include 'COMMON.GEO'
2283       include 'COMMON.LOCAL'
2284       include 'COMMON.INTERACT'
2285       include 'COMMON.IOUNITS'
2286       include 'COMMON.NAMES'
2287       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2288       common /stochcalc/ stochforcvec
2289 c
2290 c Compute the stochastic forces which contribute to velocity change
2291 c
2292       call stochastic_force(stochforcvecV)
2293 #ifndef LANG0
2294       do i=1,dimen
2295         ddt1=0.0d0
2296         ddt2=0.0d0
2297         do j=1,dimen
2298
2299           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2300 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2301           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2302         enddo
2303         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2304       enddo
2305 #endif
2306       do j=1,3
2307         d_t(j,0)=d_t_work(j)
2308       enddo
2309       ind=3
2310       do i=nnt,nct-1
2311         do j=1,3
2312           d_t(j,i)=d_t_work(ind+j)
2313         enddo
2314         ind=ind+3
2315       enddo
2316       do i=nnt,nct
2317         if (itype(i).ne.10) then
2318           inres=i+nres
2319           do j=1,3
2320             d_t(j,inres)=d_t_work(ind+j)
2321           enddo
2322           ind=ind+3
2323         endif
2324       enddo
2325       return
2326       end
2327 #endif