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

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