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(itype(i))*dsqrt(gamsc(itype(i)))
1437         enddo
1438       endif
1439 c Open the pdb file for snapshotshots
1440 #ifdef MPI
1441       if(mdpdb) then
1442         if (ilen(tmpdir).gt.0)
1443      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1444      &      liczba(:ilen(liczba))//".pdb")
1445         open(ipdb,
1446      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1447      &  //".pdb")
1448       else
1449 #ifdef NOXDR
1450         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1451      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1452      &      liczba(:ilen(liczba))//".x")
1453         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1454      &  //".x"
1455 #else
1456         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1457      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1458      &      liczba(:ilen(liczba))//".cx")
1459         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1460      &  //".cx"
1461 #endif
1462       endif
1463 #else
1464       if(mdpdb) then
1465          if (ilen(tmpdir).gt.0)
1466      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1467          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1468       else
1469          if (ilen(tmpdir).gt.0)
1470      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1471          cartname=prefix(:ilen(prefix))//"_MD.cx"
1472       endif
1473 #endif
1474       if (usampl) then
1475         write (qstr,'(256(1h ))')
1476         ipos=1
1477         do i=1,nfrag
1478           iq = qinfrag(i,iset)*10
1479           iw = wfrag(i,iset)/100
1480           if (iw.gt.0) then
1481             if(me.eq.king.or..not.out1file)
1482      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1483             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1484             ipos=ipos+7
1485           endif
1486         enddo
1487         do i=1,npair
1488           iq = qinpair(i,iset)*10
1489           iw = wpair(i,iset)/100
1490           if (iw.gt.0) then
1491             if(me.eq.king.or..not.out1file)
1492      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1493             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1494             ipos=ipos+7
1495           endif
1496         enddo
1497 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1498 #ifdef NOXDR
1499 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1500 #else
1501 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1502 #endif
1503 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1504       endif
1505       icg=1
1506       if (rest) then
1507        if (restart1file) then
1508          if (me.eq.king)
1509      &     inquire(file=mremd_rst_name,exist=file_exist)
1510            write (*,*) me," Before broadcast: file_exist",file_exist
1511          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1512      &          IERR)
1513          write (*,*) me," After broadcast: file_exist",file_exist
1514 c        inquire(file=mremd_rst_name,exist=file_exist)
1515         if(me.eq.king.or..not.out1file)
1516      &   write(iout,*) "Initial state read by master and distributed"
1517        else
1518          if (ilen(tmpdir).gt.0)
1519      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1520      &      //liczba(:ilen(liczba))//'.rst')
1521         inquire(file=rest2name,exist=file_exist)
1522        endif
1523        if(file_exist) then
1524          if(.not.restart1file) then
1525            if(me.eq.king.or..not.out1file)
1526      &      write(iout,*) "Initial state will be read from file ",
1527      &      rest2name(:ilen(rest2name))
1528            call readrst
1529          endif
1530          call rescale_weights(t_bath)
1531        else
1532         if(me.eq.king.or..not.out1file)then
1533          if (restart1file) then
1534           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1535      &       " does not exist"
1536          else
1537           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1538      &       " does not exist"
1539          endif
1540          write(iout,*) "Initial velocities randomly generated"
1541         endif
1542         call random_vel
1543         totT=0.0d0
1544        endif
1545       else
1546 c Generate initial velocities
1547         if(me.eq.king.or..not.out1file)
1548      &   write(iout,*) "Initial velocities randomly generated"
1549         call random_vel
1550         totT=0.0d0
1551       endif
1552 c      rest2name = prefix(:ilen(prefix))//'.rst'
1553       if(me.eq.king.or..not.out1file)then
1554        write (iout,*) "Initial velocities"
1555        do i=0,nres
1556          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1557      &   (d_t(j,i+nres),j=1,3)
1558        enddo
1559 c  Zeroing the total angular momentum of the system
1560        write(iout,*) "Calling the zero-angular
1561      & momentum subroutine"
1562       endif
1563       call inertia_tensor
1564 c  Getting the potential energy and forces and velocities and accelerations
1565       call vcm_vel(vcm)
1566 c      write (iout,*) "velocity of the center of the mass:"
1567 c      write (iout,*) (vcm(j),j=1,3)
1568       do j=1,3
1569         d_t(j,0)=d_t(j,0)-vcm(j)
1570       enddo
1571 c Removing the velocity of the center of mass
1572       call vcm_vel(vcm)
1573       if(me.eq.king.or..not.out1file)then
1574        write (iout,*) "vcm right after adjustment:"
1575        write (iout,*) (vcm(j),j=1,3)
1576       endif
1577       if (.not.rest) then
1578          call chainbuild
1579          if(iranconf.ne.0) then
1580           if (overlapsc) then
1581            print *, 'Calling OVERLAP_SC'
1582            call overlap_sc(fail)
1583           endif
1584
1585           if (searchsc) then
1586            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1587            print *,'SC_move',nft_sc,etot
1588            if(me.eq.king.or..not.out1file)
1589      &      write(iout,*) 'SC_move',nft_sc,etot
1590           endif
1591
1592           if(dccart)then
1593            print *, 'Calling MINIM_DC'
1594            call minim_dc(etot,iretcode,nfun)
1595           else
1596            call geom_to_var(nvar,varia)
1597            print *,'Calling MINIMIZE.'
1598            call minimize(etot,varia,iretcode,nfun)
1599            call var_to_geom(nvar,varia)
1600           endif
1601           if(me.eq.king.or..not.out1file)
1602      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1603          endif
1604       endif
1605       call chainbuild_cart
1606       call kinetic(EK)
1607       if (tbf) then
1608         call verlet_bath(EK)
1609       endif
1610       kinetic_T=2.0d0/(dimen3*Rb)*EK
1611       if(me.eq.king.or..not.out1file)then
1612        call cartprint
1613        call intout
1614       endif
1615 #ifdef MPI
1616       tt0=MPI_Wtime()
1617 #else
1618       tt0=tcpu()
1619 #endif
1620       call zerograd
1621       call etotal(potEcomp)
1622       if (large) call enerprint(potEcomp)
1623 #ifdef TIMING_ENE
1624 #ifdef MPI
1625       t_etotal=t_etotal+MPI_Wtime()-tt0
1626 #else
1627       t_etotal=t_etotal+tcpu()-tt0
1628 #endif
1629 #endif
1630       potE=potEcomp(0)
1631       call cartgrad
1632       call lagrangian
1633       call max_accel
1634       if (amax*d_time .gt. dvmax) then
1635         d_time=d_time*dvmax/amax
1636         if(me.eq.king.or..not.out1file) write (iout,*)
1637      &   "Time step reduced to",d_time,
1638      &   " because of too large initial acceleration."
1639       endif
1640       if(me.eq.king.or..not.out1file)then
1641        write(iout,*) "Potential energy and its components"
1642        call enerprint(potEcomp)
1643 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1644       endif
1645       potE=potEcomp(0)-potEcomp(20)
1646       totE=EK+potE
1647       itime=0
1648       if (ntwe.ne.0) call statout(itime)
1649       if(me.eq.king.or..not.out1file)
1650      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1651      &   " Kinetic energy",EK," Potential energy",potE,
1652      &   " Total energy",totE," Maximum acceleration ",
1653      &   amax
1654       if (large) then
1655         write (iout,*) "Initial coordinates"
1656         do i=1,nres
1657           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1658      &    (c(j,i+nres),j=1,3)
1659         enddo
1660         write (iout,*) "Initial dC"
1661         do i=0,nres
1662           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1663      &    (dc(j,i+nres),j=1,3)
1664         enddo
1665         write (iout,*) "Initial velocities"
1666         write (iout,"(13x,' backbone ',23x,' side chain')")
1667         do i=0,nres
1668           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1669      &    (d_t(j,i+nres),j=1,3)
1670         enddo
1671         write (iout,*) "Initial accelerations"
1672         do i=0,nres
1673 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1674           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1675      &    (d_a(j,i+nres),j=1,3)
1676         enddo
1677       endif
1678       do i=0,2*nres
1679         do j=1,3
1680           dc_old(j,i)=dc(j,i)
1681           d_t_old(j,i)=d_t(j,i)
1682           d_a_old(j,i)=d_a(j,i)
1683         enddo
1684 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1685       enddo
1686       if (RESPA) then
1687 #ifdef MPI
1688         tt0 =MPI_Wtime()
1689 #else
1690         tt0 = tcpu()
1691 #endif
1692         call zerograd
1693         call etotal_short(energia_short)
1694         if (large) call enerprint(potEcomp)
1695 #ifdef TIMING_ENE
1696 #ifdef MPI
1697         t_eshort=t_eshort+MPI_Wtime()-tt0
1698 #else
1699         t_eshort=t_eshort+tcpu()-tt0
1700 #endif
1701 #endif
1702         call cartgrad
1703         call lagrangian
1704         if(.not.out1file .and. large) then
1705           write (iout,*) "energia_long",energia_long(0),
1706      &     " energia_short",energia_short(0),
1707      &     " total",energia_long(0)+energia_short(0)
1708           write (iout,*) "Initial fast-force accelerations"
1709           do i=0,nres
1710             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1711      &      (d_a(j,i+nres),j=1,3)
1712           enddo
1713         endif
1714 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1715         do i=0,2*nres
1716           do j=1,3
1717             d_a_short(j,i)=d_a(j,i)
1718           enddo
1719         enddo
1720 #ifdef MPI
1721         tt0=MPI_Wtime()
1722 #else
1723         tt0=tcpu()
1724 #endif
1725         call zerograd
1726         call etotal_long(energia_long)
1727         if (large) call enerprint(potEcomp)
1728 #ifdef TIMING_ENE
1729 #ifdef MPI
1730         t_elong=t_elong+MPI_Wtime()-tt0
1731 #else
1732         t_elong=t_elong+tcpu()-tt0
1733 #endif
1734 #endif
1735         call cartgrad
1736         call lagrangian
1737         if(.not.out1file .and. large) then
1738           write (iout,*) "energia_long",energia_long(0)
1739           write (iout,*) "Initial slow-force accelerations"
1740           do i=0,nres
1741             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1742      &      (d_a(j,i+nres),j=1,3)
1743           enddo
1744         endif
1745 #ifdef MPI
1746         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1747 #else
1748         t_enegrad=t_enegrad+tcpu()-tt0
1749 #endif
1750       endif
1751       return
1752       end
1753 c-----------------------------------------------------------
1754       subroutine random_vel
1755       implicit real*8 (a-h,o-z)
1756       include 'DIMENSIONS'
1757       include 'COMMON.CONTROL'
1758       include 'COMMON.VAR'
1759       include 'COMMON.MD'
1760 #ifndef LANG0
1761       include 'COMMON.LANGEVIN'
1762 #else
1763       include 'COMMON.LANGEVIN.lang0'
1764 #endif
1765       include 'COMMON.CHAIN'
1766       include 'COMMON.DERIV'
1767       include 'COMMON.GEO'
1768       include 'COMMON.LOCAL'
1769       include 'COMMON.INTERACT'
1770       include 'COMMON.IOUNITS'
1771       include 'COMMON.NAMES'
1772       include 'COMMON.TIME1'
1773       double precision xv,sigv,lowb,highb
1774 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1775 c First generate velocities in the eigenspace of the G matrix
1776 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1777 c      call flush(iout)
1778       xv=0.0d0
1779       ii=0
1780       do i=1,dimen
1781         do k=1,3
1782           ii=ii+1
1783           sigv=dsqrt((Rb*t_bath)/geigen(i))
1784           lowb=-5*sigv
1785           highb=5*sigv
1786           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1787 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1788 c     &      " d_t_work_new",d_t_work_new(ii)
1789         enddo
1790       enddo
1791 c diagnostics
1792 c      Ek1=0.0d0
1793 c      ii=0
1794 c      do i=1,dimen
1795 c        do k=1,3
1796 c          ii=ii+1
1797 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1798 c        enddo
1799 c      enddo
1800 c      write (iout,*) "Ek from eigenvectors",Ek1
1801 c end diagnostics
1802 c Transform velocities to UNRES coordinate space
1803       do k=0,2
1804         do i=1,dimen
1805           ind=(i-1)*3+k+1
1806           d_t_work(ind)=0.0d0
1807           do j=1,dimen
1808             d_t_work(ind)=d_t_work(ind)
1809      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1810           enddo
1811 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1812 c          call flush(iout)
1813         enddo
1814       enddo
1815 c Transfer to the d_t vector
1816       do j=1,3
1817         d_t(j,0)=d_t_work(j)
1818       enddo
1819       ind=3
1820       do i=nnt,nct-1
1821         do j=1,3
1822           ind=ind+1
1823           d_t(j,i)=d_t_work(ind)
1824         enddo
1825       enddo
1826       do i=nnt,nct
1827         if (itype(i).ne.10 .and. itype(i).ne.21) then
1828           do j=1,3
1829             ind=ind+1
1830             d_t(j,i+nres)=d_t_work(ind)
1831           enddo
1832         endif
1833       enddo
1834 c      call kinetic(EK)
1835 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1836 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1837 c      call flush(iout)
1838       return
1839       end
1840 #ifndef LANG0
1841 c-----------------------------------------------------------
1842       subroutine sd_verlet_p_setup
1843 c Sets up the parameters of stochastic Verlet algorithm
1844       implicit real*8 (a-h,o-z)
1845       include 'DIMENSIONS'
1846 #ifdef MPI
1847       include 'mpif.h'
1848 #endif
1849       include 'COMMON.CONTROL'
1850       include 'COMMON.VAR'
1851       include 'COMMON.MD'
1852 #ifndef LANG0
1853       include 'COMMON.LANGEVIN'
1854 #else
1855       include 'COMMON.LANGEVIN.lang0'
1856 #endif
1857       include 'COMMON.CHAIN'
1858       include 'COMMON.DERIV'
1859       include 'COMMON.GEO'
1860       include 'COMMON.LOCAL'
1861       include 'COMMON.INTERACT'
1862       include 'COMMON.IOUNITS'
1863       include 'COMMON.NAMES'
1864       include 'COMMON.TIME1'
1865       double precision emgdt(MAXRES6),
1866      & pterm,vterm,rho,rhoc,vsig,
1867      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1868      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1869      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1870       logical lprn /.false./
1871       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1872       double precision ktm
1873 #ifdef MPI
1874       tt0 = MPI_Wtime()
1875 #else
1876       tt0 = tcpu()
1877 #endif
1878 c
1879 c AL 8/17/04 Code adapted from tinker
1880 c
1881 c Get the frictional and random terms for stochastic dynamics in the
1882 c eigenspace of mass-scaled UNRES friction matrix
1883 c
1884       do i = 1, dimen
1885             gdt = fricgam(i) * d_time
1886 c
1887 c Stochastic dynamics reduces to simple MD for zero friction
1888 c
1889             if (gdt .le. zero) then
1890                pfric_vec(i) = 1.0d0
1891                vfric_vec(i) = d_time
1892                afric_vec(i) = 0.5d0 * d_time * d_time
1893                prand_vec(i) = 0.0d0
1894                vrand_vec1(i) = 0.0d0
1895                vrand_vec2(i) = 0.0d0
1896 c
1897 c Analytical expressions when friction coefficient is large
1898 c
1899             else
1900                if (gdt .ge. gdt_radius) then
1901                   egdt = dexp(-gdt)
1902                   pfric_vec(i) = egdt
1903                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1904                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1905                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1906                   vterm = 1.0d0 - egdt**2
1907                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1908 c
1909 c Use series expansions when friction coefficient is small
1910 c
1911                else
1912                   gdt2 = gdt * gdt
1913                   gdt3 = gdt * gdt2
1914                   gdt4 = gdt2 * gdt2
1915                   gdt5 = gdt2 * gdt3
1916                   gdt6 = gdt3 * gdt3
1917                   gdt7 = gdt3 * gdt4
1918                   gdt8 = gdt4 * gdt4
1919                   gdt9 = gdt4 * gdt5
1920                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1921      &                          - gdt5/120.0d0 + gdt6/720.0d0
1922      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
1923      &                          - gdt9/362880.0d0) / fricgam(i)**2
1924                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1925                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1926                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1927      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1928      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1929      &                       + 127.0d0*gdt9/90720.0d0
1930                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1931      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1932      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1933      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1934                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1935      &                       - 17.0d0*gdt2/1280.0d0
1936      &                       + 17.0d0*gdt3/6144.0d0
1937      &                       + 40967.0d0*gdt4/34406400.0d0
1938      &                       - 57203.0d0*gdt5/275251200.0d0
1939      &                       - 1429487.0d0*gdt6/13212057600.0d0)
1940                end if
1941 c
1942 c Compute the scaling factors of random terms for the nonzero friction case
1943 c
1944                ktm = 0.5d0*d_time/fricgam(i)
1945                psig = dsqrt(ktm*pterm) / fricgam(i)
1946                vsig = dsqrt(ktm*vterm)
1947                rhoc = dsqrt(1.0d0 - rho*rho)
1948                prand_vec(i) = psig
1949                vrand_vec1(i) = vsig * rho
1950                vrand_vec2(i) = vsig * rhoc
1951             end if
1952       end do
1953       if (lprn) then
1954       write (iout,*)
1955      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1956      &  " vrand_vec2"
1957       do i=1,dimen
1958         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1959      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1960       enddo
1961       endif
1962 c
1963 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1964 c
1965 #ifndef   LANG0
1966       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1967       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1968       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1969       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1970       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1971       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1972 #endif
1973 #ifdef MPI
1974       t_sdsetup=t_sdsetup+MPI_Wtime()
1975 #else
1976       t_sdsetup=t_sdsetup+tcpu()-tt0
1977 #endif
1978       return
1979       end
1980 c-------------------------------------------------------------
1981       subroutine eigtransf1(n,ndim,ab,d,c)
1982       implicit none
1983       integer n,ndim
1984       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1985       integer i,j,k
1986       do i=1,n
1987         do j=1,n
1988           c(i,j)=0.0d0
1989           do k=1,n
1990             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1991           enddo
1992         enddo
1993       enddo
1994       return
1995       end
1996 c-------------------------------------------------------------
1997       subroutine eigtransf(n,ndim,a,b,d,c)
1998       implicit none
1999       integer n,ndim
2000       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2001       integer i,j,k
2002       do i=1,n
2003         do j=1,n
2004           c(i,j)=0.0d0
2005           do k=1,n
2006             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2007           enddo
2008         enddo
2009       enddo
2010       return
2011       end
2012 c-------------------------------------------------------------
2013       subroutine sd_verlet1
2014 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2015       implicit real*8 (a-h,o-z)
2016       include 'DIMENSIONS'
2017       include 'COMMON.CONTROL'
2018       include 'COMMON.VAR'
2019       include 'COMMON.MD'
2020 #ifndef LANG0
2021       include 'COMMON.LANGEVIN'
2022 #else
2023       include 'COMMON.LANGEVIN.lang0'
2024 #endif
2025       include 'COMMON.CHAIN'
2026       include 'COMMON.DERIV'
2027       include 'COMMON.GEO'
2028       include 'COMMON.LOCAL'
2029       include 'COMMON.INTERACT'
2030       include 'COMMON.IOUNITS'
2031       include 'COMMON.NAMES'
2032       double precision stochforcvec(MAXRES6)
2033       common /stochcalc/ stochforcvec
2034       logical lprn /.false./
2035
2036 c      write (iout,*) "dc_old"
2037 c      do i=0,nres
2038 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2039 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2040 c      enddo
2041       do j=1,3
2042         dc_work(j)=dc_old(j,0)
2043         d_t_work(j)=d_t_old(j,0)
2044         d_a_work(j)=d_a_old(j,0)
2045       enddo
2046       ind=3
2047       do i=nnt,nct-1
2048         do j=1,3
2049           dc_work(ind+j)=dc_old(j,i)
2050           d_t_work(ind+j)=d_t_old(j,i)
2051           d_a_work(ind+j)=d_a_old(j,i)
2052         enddo
2053         ind=ind+3
2054       enddo
2055       do i=nnt,nct
2056         if (itype(i).ne.10 .and. itype(i).ne.21) then
2057           do j=1,3
2058             dc_work(ind+j)=dc_old(j,i+nres)
2059             d_t_work(ind+j)=d_t_old(j,i+nres)
2060             d_a_work(ind+j)=d_a_old(j,i+nres)
2061           enddo
2062           ind=ind+3
2063         endif
2064       enddo
2065 #ifndef LANG0
2066       if (lprn) then
2067       write (iout,*)
2068      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2069      &  " vrand_mat2"
2070       do i=1,dimen
2071         do j=1,dimen
2072           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2073      &      vfric_mat(i,j),afric_mat(i,j),
2074      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2075         enddo
2076       enddo
2077       endif
2078       do i=1,dimen
2079         ddt1=0.0d0
2080         ddt2=0.0d0
2081         do j=1,dimen
2082           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2083      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2084           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2085           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2086         enddo
2087         d_t_work_new(i)=ddt1+0.5d0*ddt2
2088         d_t_work(i)=ddt1+ddt2
2089       enddo
2090 #endif
2091       do j=1,3
2092         dc(j,0)=dc_work(j)
2093         d_t(j,0)=d_t_work(j)
2094       enddo
2095       ind=3
2096       do i=nnt,nct-1
2097         do j=1,3
2098           dc(j,i)=dc_work(ind+j)
2099           d_t(j,i)=d_t_work(ind+j)
2100         enddo
2101         ind=ind+3
2102       enddo
2103       do i=nnt,nct
2104         if (itype(i).ne.10) then
2105           inres=i+nres
2106           do j=1,3
2107             dc(j,inres)=dc_work(ind+j)
2108             d_t(j,inres)=d_t_work(ind+j)
2109           enddo
2110           ind=ind+3
2111         endif
2112       enddo
2113       return
2114       end
2115 c--------------------------------------------------------------------------
2116       subroutine sd_verlet2
2117 c  Calculating the adjusted velocities for accelerations
2118       implicit real*8 (a-h,o-z)
2119       include 'DIMENSIONS'
2120       include 'COMMON.CONTROL'
2121       include 'COMMON.VAR'
2122       include 'COMMON.MD'
2123 #ifndef LANG0
2124       include 'COMMON.LANGEVIN'
2125 #else
2126       include 'COMMON.LANGEVIN.lang0'
2127 #endif
2128       include 'COMMON.CHAIN'
2129       include 'COMMON.DERIV'
2130       include 'COMMON.GEO'
2131       include 'COMMON.LOCAL'
2132       include 'COMMON.INTERACT'
2133       include 'COMMON.IOUNITS'
2134       include 'COMMON.NAMES'
2135       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2136       common /stochcalc/ stochforcvec
2137 c
2138 c Compute the stochastic forces which contribute to velocity change
2139 c
2140       call stochastic_force(stochforcvecV)
2141
2142 #ifndef LANG0
2143       do i=1,dimen
2144         ddt1=0.0d0
2145         ddt2=0.0d0
2146         do j=1,dimen
2147           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2148           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2149      &     vrand_mat2(i,j)*stochforcvecV(j)
2150         enddo
2151         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2152       enddo
2153 #endif
2154       do j=1,3
2155         d_t(j,0)=d_t_work(j)
2156       enddo
2157       ind=3
2158       do i=nnt,nct-1
2159         do j=1,3
2160           d_t(j,i)=d_t_work(ind+j)
2161         enddo
2162         ind=ind+3
2163       enddo
2164       do i=nnt,nct
2165         if (itype(i).ne.10 .and. itype(i).ne.21) then
2166           inres=i+nres
2167           do j=1,3
2168             d_t(j,inres)=d_t_work(ind+j)
2169           enddo
2170           ind=ind+3
2171         endif
2172       enddo
2173       return
2174       end
2175 c-----------------------------------------------------------
2176       subroutine sd_verlet_ciccotti_setup
2177 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2178 c version
2179       implicit real*8 (a-h,o-z)
2180       include 'DIMENSIONS'
2181 #ifdef MPI
2182       include 'mpif.h'
2183 #endif
2184       include 'COMMON.CONTROL'
2185       include 'COMMON.VAR'
2186       include 'COMMON.MD'
2187 #ifndef LANG0
2188       include 'COMMON.LANGEVIN'
2189 #else
2190       include 'COMMON.LANGEVIN.lang0'
2191 #endif
2192       include 'COMMON.CHAIN'
2193       include 'COMMON.DERIV'
2194       include 'COMMON.GEO'
2195       include 'COMMON.LOCAL'
2196       include 'COMMON.INTERACT'
2197       include 'COMMON.IOUNITS'
2198       include 'COMMON.NAMES'
2199       include 'COMMON.TIME1'
2200       double precision emgdt(MAXRES6),
2201      & pterm,vterm,rho,rhoc,vsig,
2202      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2203      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2204      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2205       logical lprn /.false./
2206       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2207       double precision ktm
2208 #ifdef MPI
2209       tt0 = MPI_Wtime()
2210 #else
2211       tt0 = tcpu()
2212 #endif
2213 c
2214 c AL 8/17/04 Code adapted from tinker
2215 c
2216 c Get the frictional and random terms for stochastic dynamics in the
2217 c eigenspace of mass-scaled UNRES friction matrix
2218 c
2219       do i = 1, dimen
2220             write (iout,*) "i",i," fricgam",fricgam(i)
2221             gdt = fricgam(i) * d_time
2222 c
2223 c Stochastic dynamics reduces to simple MD for zero friction
2224 c
2225             if (gdt .le. zero) then
2226                pfric_vec(i) = 1.0d0
2227                vfric_vec(i) = d_time
2228                afric_vec(i) = 0.5d0*d_time*d_time
2229                prand_vec(i) = afric_vec(i)
2230                vrand_vec2(i) = vfric_vec(i)
2231 c
2232 c Analytical expressions when friction coefficient is large
2233 c
2234             else
2235                egdt = dexp(-gdt)
2236                pfric_vec(i) = egdt
2237                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2238                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2239                prand_vec(i) = afric_vec(i)
2240                vrand_vec2(i) = vfric_vec(i)
2241 c
2242 c Compute the scaling factors of random terms for the nonzero friction case
2243 c
2244 c               ktm = 0.5d0*d_time/fricgam(i)
2245 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2246 c               vsig = dsqrt(ktm*vterm)
2247 c               prand_vec(i) = psig*afric_vec(i)
2248 c               vrand_vec2(i) = vsig*vfric_vec(i)
2249             end if
2250       end do
2251       if (lprn) then
2252       write (iout,*)
2253      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2254      &  " vrand_vec2"
2255       do i=1,dimen
2256         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2257      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2258       enddo
2259       endif
2260 c
2261 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2262 c
2263       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2264       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2265       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2266       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2267       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2268 #ifdef MPI
2269       t_sdsetup=t_sdsetup+MPI_Wtime()
2270 #else
2271       t_sdsetup=t_sdsetup+tcpu()-tt0
2272 #endif
2273       return
2274       end
2275 c-------------------------------------------------------------
2276       subroutine sd_verlet1_ciccotti
2277 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2278       implicit real*8 (a-h,o-z)
2279       include 'DIMENSIONS'
2280 #ifdef MPI
2281       include 'mpif.h'
2282 #endif
2283       include 'COMMON.CONTROL'
2284       include 'COMMON.VAR'
2285       include 'COMMON.MD'
2286 #ifndef LANG0
2287       include 'COMMON.LANGEVIN'
2288 #else
2289       include 'COMMON.LANGEVIN.lang0'
2290 #endif
2291       include 'COMMON.CHAIN'
2292       include 'COMMON.DERIV'
2293       include 'COMMON.GEO'
2294       include 'COMMON.LOCAL'
2295       include 'COMMON.INTERACT'
2296       include 'COMMON.IOUNITS'
2297       include 'COMMON.NAMES'
2298       double precision stochforcvec(MAXRES6)
2299       common /stochcalc/ stochforcvec
2300       logical lprn /.false./
2301
2302 c      write (iout,*) "dc_old"
2303 c      do i=0,nres
2304 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2305 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2306 c      enddo
2307       do j=1,3
2308         dc_work(j)=dc_old(j,0)
2309         d_t_work(j)=d_t_old(j,0)
2310         d_a_work(j)=d_a_old(j,0)
2311       enddo
2312       ind=3
2313       do i=nnt,nct-1
2314         do j=1,3
2315           dc_work(ind+j)=dc_old(j,i)
2316           d_t_work(ind+j)=d_t_old(j,i)
2317           d_a_work(ind+j)=d_a_old(j,i)
2318         enddo
2319         ind=ind+3
2320       enddo
2321       do i=nnt,nct
2322         if (itype(i).ne.10 .and. itype(i).ne.21) then
2323           do j=1,3
2324             dc_work(ind+j)=dc_old(j,i+nres)
2325             d_t_work(ind+j)=d_t_old(j,i+nres)
2326             d_a_work(ind+j)=d_a_old(j,i+nres)
2327           enddo
2328           ind=ind+3
2329         endif
2330       enddo
2331
2332 #ifndef LANG0
2333       if (lprn) then
2334       write (iout,*)
2335      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2336      &  " vrand_mat2"
2337       do i=1,dimen
2338         do j=1,dimen
2339           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2340      &      vfric_mat(i,j),afric_mat(i,j),
2341      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2342         enddo
2343       enddo
2344       endif
2345       do i=1,dimen
2346         ddt1=0.0d0
2347         ddt2=0.0d0
2348         do j=1,dimen
2349           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2350      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2351           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2352           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2353         enddo
2354         d_t_work_new(i)=ddt1+0.5d0*ddt2
2355         d_t_work(i)=ddt1+ddt2
2356       enddo
2357 #endif
2358       do j=1,3
2359         dc(j,0)=dc_work(j)
2360         d_t(j,0)=d_t_work(j)
2361       enddo
2362       ind=3
2363       do i=nnt,nct-1
2364         do j=1,3
2365           dc(j,i)=dc_work(ind+j)
2366           d_t(j,i)=d_t_work(ind+j)
2367         enddo
2368         ind=ind+3
2369       enddo
2370       do i=nnt,nct
2371         if (itype(i).ne.10 .and. itype(i).ne.21) then
2372           inres=i+nres
2373           do j=1,3
2374             dc(j,inres)=dc_work(ind+j)
2375             d_t(j,inres)=d_t_work(ind+j)
2376           enddo
2377           ind=ind+3
2378         endif
2379       enddo
2380       return
2381       end
2382 c--------------------------------------------------------------------------
2383       subroutine sd_verlet2_ciccotti
2384 c  Calculating the adjusted velocities for accelerations
2385       implicit real*8 (a-h,o-z)
2386       include 'DIMENSIONS'
2387       include 'COMMON.CONTROL'
2388       include 'COMMON.VAR'
2389       include 'COMMON.MD'
2390 #ifndef LANG0
2391       include 'COMMON.LANGEVIN'
2392 #else
2393       include 'COMMON.LANGEVIN.lang0'
2394 #endif
2395       include 'COMMON.CHAIN'
2396       include 'COMMON.DERIV'
2397       include 'COMMON.GEO'
2398       include 'COMMON.LOCAL'
2399       include 'COMMON.INTERACT'
2400       include 'COMMON.IOUNITS'
2401       include 'COMMON.NAMES'
2402       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2403       common /stochcalc/ stochforcvec
2404 c
2405 c Compute the stochastic forces which contribute to velocity change
2406 c
2407       call stochastic_force(stochforcvecV)
2408 #ifndef LANG0
2409       do i=1,dimen
2410         ddt1=0.0d0
2411         ddt2=0.0d0
2412         do j=1,dimen
2413
2414           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2415 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2416           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2417         enddo
2418         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2419       enddo
2420 #endif
2421       do j=1,3
2422         d_t(j,0)=d_t_work(j)
2423       enddo
2424       ind=3
2425       do i=nnt,nct-1
2426         do j=1,3
2427           d_t(j,i)=d_t_work(ind+j)
2428         enddo
2429         ind=ind+3
2430       enddo
2431       do i=nnt,nct
2432         if (itype(i).ne.10 .and. itype(i).ne.21) then
2433           inres=i+nres
2434           do j=1,3
2435             d_t(j,inres)=d_t_work(ind+j)
2436           enddo
2437           ind=ind+3
2438         endif
2439       enddo
2440       return
2441       end
2442 #endif