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