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