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 #ifdef MPI
1091       time00=MPI_Wtime()
1092 #else
1093       time00=tcpu()
1094 #endif
1095       call friction_force
1096 #ifdef MPI
1097       time_fric=time_fric+MPI_Wtime()-time00
1098       time00=MPI_Wtime()
1099 #else
1100       time_fric=time_fric+tcpu()-time00
1101       time00=tcpu()
1102 #endif
1103       call stochastic_force(stochforcvec)
1104 #ifdef MPI
1105       time_stoch=time_stoch+MPI_Wtime()-time00
1106 #else
1107       time_stoch=time_stoch+tcpu()-time00
1108 #endif
1109 c
1110 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1111 c forces (d_as_work)
1112 c
1113       call ginv_mult(fric_work, d_af_work)
1114       call ginv_mult(stochforcvec, d_as_work)
1115       return
1116       end
1117 c---------------------------------------------------------------------
1118       subroutine sddir_verlet1
1119 c Applying velocity Verlet algorithm - step 1 to velocities
1120       implicit real*8 (a-h,o-z)
1121       include 'DIMENSIONS'
1122       include 'COMMON.CONTROL'
1123       include 'COMMON.VAR'
1124       include 'COMMON.MD'
1125 #ifndef LANG0
1126       include 'COMMON.LANGEVIN'
1127 #else
1128       include 'COMMON.LANGEVIN.lang0'
1129 #endif
1130       include 'COMMON.CHAIN'
1131       include 'COMMON.DERIV'
1132       include 'COMMON.GEO'
1133       include 'COMMON.LOCAL'
1134       include 'COMMON.INTERACT'
1135       include 'COMMON.IOUNITS'
1136       include 'COMMON.NAMES'
1137 c Revised 3/31/05 AL: correlation between random contributions to
1138 c position and velocity increments included.
1139       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1140       double precision adt,adt2
1141 c
1142 c Add the contribution from BOTH friction and stochastic force to the
1143 c coordinates, but ONLY the contribution from the friction forces to velocities
1144 c
1145       do j=1,3
1146         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1147         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1148         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1149         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1150         d_t(j,0)=d_t_old(j,0)+adt
1151       enddo
1152       ind=3
1153       do i=nnt,nct-1
1154         do j=1,3
1155           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1156           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1157           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1158           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1159           d_t(j,i)=d_t_old(j,i)+adt
1160         enddo
1161         ind=ind+3
1162       enddo
1163       do i=nnt,nct
1164         if (itype(i).ne.10 .and. itype(i).ne.21) then
1165           inres=i+nres
1166           do j=1,3
1167             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1168             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1169             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1170             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1171             d_t(j,inres)=d_t_old(j,inres)+adt
1172           enddo
1173           ind=ind+3
1174         endif
1175       enddo
1176       return
1177       end
1178 c---------------------------------------------------------------------
1179       subroutine sddir_verlet2
1180 c  Calculating the adjusted velocities for accelerations
1181       implicit real*8 (a-h,o-z)
1182       include 'DIMENSIONS'
1183       include 'COMMON.CONTROL'
1184       include 'COMMON.VAR'
1185       include 'COMMON.MD'
1186 #ifndef LANG0
1187       include 'COMMON.LANGEVIN'
1188 #else
1189       include 'COMMON.LANGEVIN.lang0'
1190 #endif
1191       include 'COMMON.CHAIN'
1192       include 'COMMON.DERIV'
1193       include 'COMMON.GEO'
1194       include 'COMMON.LOCAL'
1195       include 'COMMON.INTERACT'
1196       include 'COMMON.IOUNITS'
1197       include 'COMMON.NAMES'
1198       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1199       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1200 c Revised 3/31/05 AL: correlation between random contributions to
1201 c position and velocity increments included.
1202 c The correlation coefficients are calculated at low-friction limit.
1203 c Also, friction forces are now not calculated with new velocities.
1204
1205 c      call friction_force
1206       call stochastic_force(stochforcvec)
1207 c
1208 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1209 c forces (d_as_work)
1210 c
1211       call ginv_mult(stochforcvec, d_as_work1)
1212
1213 c
1214 c Update velocities
1215 c
1216       do j=1,3
1217         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1218      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1219       enddo
1220       ind=3
1221       do i=nnt,nct-1
1222         do j=1,3
1223           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1224      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1225         enddo
1226         ind=ind+3
1227       enddo
1228       do i=nnt,nct
1229         if (itype(i).ne.10 .and. itype(i).ne.21) then
1230           inres=i+nres
1231           do j=1,3
1232             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1233      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1234      &       +cos60*d_as_work1(ind+j))*d_time
1235           enddo
1236           ind=ind+3
1237         endif
1238       enddo
1239       return
1240       end
1241 c---------------------------------------------------------------------
1242       subroutine max_accel
1243 c
1244 c Find the maximum difference in the accelerations of the the sites
1245 c at the beginning and the end of the time step.
1246 c
1247       implicit real*8 (a-h,o-z)
1248       include 'DIMENSIONS'
1249       include 'COMMON.CONTROL'
1250       include 'COMMON.VAR'
1251       include 'COMMON.MD'
1252       include 'COMMON.CHAIN'
1253       include 'COMMON.DERIV'
1254       include 'COMMON.GEO'
1255       include 'COMMON.LOCAL'
1256       include 'COMMON.INTERACT'
1257       include 'COMMON.IOUNITS'
1258       double precision aux(3),accel(3),accel_old(3),dacc
1259       do j=1,3
1260 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1261          accel_old(j)=d_a_old(j,0)
1262          accel(j)=d_a(j,0)
1263       enddo
1264       amax=0.0d0
1265       do i=nnt,nct
1266 c Backbone
1267         if (i.lt.nct) then
1268 c 7/3/08 changed to asymmetric difference
1269           do j=1,3
1270 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1271             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1272             accel(j)=accel(j)+0.5d0*d_a(j,i)
1273 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1274             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1275               dacc=dabs(accel(j)-accel_old(j))
1276               if (dacc.gt.amax) amax=dacc
1277             endif
1278           enddo
1279         endif
1280       enddo
1281 c Side chains
1282       do j=1,3
1283 c        accel(j)=aux(j)
1284         accel_old(j)=d_a_old(j,0)
1285         accel(j)=d_a(j,0)
1286       enddo
1287       if (nnt.eq.2) then
1288         do j=1,3
1289           accel_old(j)=accel_old(j)+d_a_old(j,1)
1290           accel(j)=accel(j)+d_a(j,1)
1291         enddo
1292       endif
1293       do i=nnt,nct
1294         if (itype(i).ne.10 .and. itype(i).ne.21) then
1295           do j=1,3
1296 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1297             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1298             accel(j)=accel(j)+d_a(j,i+nres)
1299           enddo
1300         endif
1301         do j=1,3
1302 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1303           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1304             dacc=dabs(accel(j)-accel_old(j))
1305             if (dacc.gt.amax) amax=dacc
1306           endif
1307         enddo
1308         do j=1,3
1309           accel_old(j)=accel_old(j)+d_a_old(j,i)
1310           accel(j)=accel(j)+d_a(j,i)
1311 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1312         enddo
1313       enddo
1314       return
1315       end
1316 c---------------------------------------------------------------------
1317       subroutine predict_edrift(epdrift)
1318 c
1319 c Predict the drift of the potential energy
1320 c
1321       implicit real*8 (a-h,o-z)
1322       include 'DIMENSIONS'
1323       include 'COMMON.CONTROL'
1324       include 'COMMON.VAR'
1325       include 'COMMON.MD'
1326       include 'COMMON.CHAIN'
1327       include 'COMMON.DERIV'
1328       include 'COMMON.GEO'
1329       include 'COMMON.LOCAL'
1330       include 'COMMON.INTERACT'
1331       include 'COMMON.IOUNITS'
1332       include 'COMMON.MUCA'
1333       double precision epdrift,epdriftij
1334 c Drift of the potential energy
1335       epdrift=0.0d0
1336       do i=nnt,nct
1337 c Backbone
1338         if (i.lt.nct) then
1339           do j=1,3
1340             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1341             if (lmuca) epdriftij=epdriftij*factor
1342 c            write (iout,*) "back",i,j,epdriftij
1343             if (epdriftij.gt.epdrift) epdrift=epdriftij
1344           enddo
1345         endif
1346 c Side chains
1347         if (itype(i).ne.10 .and. itype(i).ne.21) then
1348           do j=1,3
1349             epdriftij=
1350      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1351             if (lmuca) epdriftij=epdriftij*factor
1352 c            write (iout,*) "side",i,j,epdriftij
1353             if (epdriftij.gt.epdrift) epdrift=epdriftij
1354           enddo
1355         endif
1356       enddo
1357       epdrift=0.5d0*epdrift*d_time*d_time
1358 c      write (iout,*) "epdrift",epdrift
1359       return
1360       end
1361 c-----------------------------------------------------------------------
1362       subroutine verlet_bath
1363 c
1364 c  Coupling to the thermostat by using the Berendsen algorithm
1365 c
1366       implicit real*8 (a-h,o-z)
1367       include 'DIMENSIONS'
1368       include 'COMMON.CONTROL'
1369       include 'COMMON.VAR'
1370       include 'COMMON.MD'
1371       include 'COMMON.CHAIN'
1372       include 'COMMON.DERIV'
1373       include 'COMMON.GEO'
1374       include 'COMMON.LOCAL'
1375       include 'COMMON.INTERACT'
1376       include 'COMMON.IOUNITS'
1377       include 'COMMON.NAMES'
1378       double precision T_half,fact
1379 c
1380       T_half=2.0d0/(dimen3*Rb)*EK
1381       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1382 c      write(iout,*) "T_half", T_half
1383 c      write(iout,*) "EK", EK
1384 c      write(iout,*) "fact", fact
1385       do j=1,3
1386         d_t(j,0)=fact*d_t(j,0)
1387       enddo
1388       do i=nnt,nct-1
1389         do j=1,3
1390           d_t(j,i)=fact*d_t(j,i)
1391         enddo
1392       enddo
1393       do i=nnt,nct
1394         if (itype(i).ne.10 .and. itype(i).ne.21) then
1395           inres=i+nres
1396           do j=1,3
1397             d_t(j,inres)=fact*d_t(j,inres)
1398           enddo
1399         endif
1400       enddo
1401       return
1402       end
1403 c---------------------------------------------------------
1404       subroutine init_MD
1405 c  Set up the initial conditions of a MD simulation
1406       implicit real*8 (a-h,o-z)
1407       include 'DIMENSIONS'
1408 #ifdef MP
1409       include 'mpif.h'
1410       character*16 form
1411       integer IERROR,ERRCODE
1412 #endif
1413       include 'COMMON.SETUP'
1414       include 'COMMON.CONTROL'
1415       include 'COMMON.VAR'
1416       include 'COMMON.MD'
1417 #ifndef LANG0
1418       include 'COMMON.LANGEVIN'
1419 #else
1420       include 'COMMON.LANGEVIN.lang0'
1421 #endif
1422       include 'COMMON.CHAIN'
1423       include 'COMMON.DERIV'
1424       include 'COMMON.GEO'
1425       include 'COMMON.LOCAL'
1426       include 'COMMON.INTERACT'
1427       include 'COMMON.IOUNITS'
1428       include 'COMMON.NAMES'
1429       include 'COMMON.REMD'
1430       real*8 energia_long(0:n_ene),
1431      &  energia_short(0:n_ene),vcm(3),incr(3)
1432       double precision cm(3),L(3),xv,sigv,lowb,highb
1433       double precision varia(maxvar)
1434       character*256 qstr
1435       integer ilen
1436       external ilen
1437       character*50 tytul
1438       logical file_exist
1439       common /gucio/ cm
1440       d_time0=d_time
1441 c      write(iout,*) "d_time", d_time
1442 c Compute the standard deviations of stochastic forces for Langevin dynamics
1443 c if the friction coefficients do not depend on surface area
1444       if (lang.gt.0 .and. .not.surfarea) then
1445         do i=nnt,nct-1
1446           stdforcp(i)=stdfp*dsqrt(gamp)
1447         enddo
1448         do i=nnt,nct
1449           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
1450         enddo
1451       endif
1452 c Open the pdb file for snapshotshots
1453 #ifdef MPI
1454       if(mdpdb) then
1455         if (ilen(tmpdir).gt.0)
1456      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1457      &      liczba(:ilen(liczba))//".pdb")
1458         open(ipdb,
1459      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1460      &  //".pdb")
1461       else
1462 #ifdef NOXDR
1463         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1464      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1465      &      liczba(:ilen(liczba))//".x")
1466         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1467      &  //".x"
1468 #else
1469         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1470      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1471      &      liczba(:ilen(liczba))//".cx")
1472         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1473      &  //".cx"
1474 #endif
1475       endif
1476 #else
1477       if(mdpdb) then
1478          if (ilen(tmpdir).gt.0)
1479      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1480          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1481       else
1482          if (ilen(tmpdir).gt.0)
1483      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1484          cartname=prefix(:ilen(prefix))//"_MD.cx"
1485       endif
1486 #endif
1487       if (usampl) then
1488         write (qstr,'(256(1h ))')
1489         ipos=1
1490         do i=1,nfrag
1491           iq = qinfrag(i,iset)*10
1492           iw = wfrag(i,iset)/100
1493           if (iw.gt.0) then
1494             if(me.eq.king.or..not.out1file)
1495      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1496             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1497             ipos=ipos+7
1498           endif
1499         enddo
1500         do i=1,npair
1501           iq = qinpair(i,iset)*10
1502           iw = wpair(i,iset)/100
1503           if (iw.gt.0) then
1504             if(me.eq.king.or..not.out1file)
1505      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1506             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1507             ipos=ipos+7
1508           endif
1509         enddo
1510 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1511 #ifdef NOXDR
1512 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1513 #else
1514 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1515 #endif
1516 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1517       endif
1518       icg=1
1519       if (rest) then
1520        if (restart1file) then
1521          if (me.eq.king)
1522      &     inquire(file=mremd_rst_name,exist=file_exist)
1523 #ifdef MPI
1524            write (*,*) me," Before broadcast: file_exist",file_exist
1525          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1526      &          IERR)
1527          write (*,*) me," After broadcast: file_exist",file_exist
1528 c        inquire(file=mremd_rst_name,exist=file_exist)
1529 #endif
1530         if(me.eq.king.or..not.out1file)
1531      &   write(iout,*) "Initial state read by master and distributed"
1532        else
1533          if (ilen(tmpdir).gt.0)
1534      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1535      &      //liczba(:ilen(liczba))//'.rst')
1536         inquire(file=rest2name,exist=file_exist)
1537        endif
1538        if(file_exist) then
1539          if(.not.restart1file) then
1540            if(me.eq.king.or..not.out1file)
1541      &      write(iout,*) "Initial state will be read from file ",
1542      &      rest2name(:ilen(rest2name))
1543            call readrst
1544          endif
1545          call rescale_weights(t_bath)
1546        else
1547         if(me.eq.king.or..not.out1file)then
1548          if (restart1file) then
1549           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1550      &       " does not exist"
1551          else
1552           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1553      &       " does not exist"
1554          endif
1555          write(iout,*) "Initial velocities randomly generated"
1556         endif
1557         call random_vel
1558         totT=0.0d0
1559        endif
1560       else
1561 c Generate initial velocities
1562         if(me.eq.king.or..not.out1file)
1563      &   write(iout,*) "Initial velocities randomly generated"
1564         call random_vel
1565         totT=0.0d0
1566       endif
1567 c      rest2name = prefix(:ilen(prefix))//'.rst'
1568       if(me.eq.king.or..not.out1file)then
1569        write (iout,*) "Initial velocities"
1570        do i=0,nres
1571          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1572      &   (d_t(j,i+nres),j=1,3)
1573        enddo
1574 c  Zeroing the total angular momentum of the system
1575        write(iout,*) "Calling the zero-angular
1576      & momentum subroutine"
1577       endif
1578       call inertia_tensor
1579 c  Getting the potential energy and forces and velocities and accelerations
1580       call vcm_vel(vcm)
1581 c      write (iout,*) "velocity of the center of the mass:"
1582 c      write (iout,*) (vcm(j),j=1,3)
1583       do j=1,3
1584         d_t(j,0)=d_t(j,0)-vcm(j)
1585       enddo
1586 c Removing the velocity of the center of mass
1587       call vcm_vel(vcm)
1588       if(me.eq.king.or..not.out1file)then
1589        write (iout,*) "vcm right after adjustment:"
1590        write (iout,*) (vcm(j),j=1,3)
1591       endif
1592       if (.not.rest) then
1593          call chainbuild
1594          if(iranconf.ne.0) then
1595           if (overlapsc) then
1596            print *, 'Calling OVERLAP_SC'
1597            call overlap_sc(fail)
1598           endif
1599
1600           if (searchsc) then
1601            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1602            print *,'SC_move',nft_sc,etot
1603            if(me.eq.king.or..not.out1file)
1604      &      write(iout,*) 'SC_move',nft_sc,etot
1605           endif
1606
1607           if(dccart)then
1608            print *, 'Calling MINIM_DC'
1609            call minim_dc(etot,iretcode,nfun)
1610           else
1611            call geom_to_var(nvar,varia)
1612            print *,'Calling MINIMIZE.'
1613            call minimize(etot,varia,iretcode,nfun)
1614            call var_to_geom(nvar,varia)
1615           endif
1616           if(me.eq.king.or..not.out1file)
1617      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1618          endif
1619       endif
1620       call chainbuild_cart
1621       call kinetic(EK)
1622       if (tbf) then
1623         call verlet_bath
1624       endif
1625       kinetic_T=2.0d0/(dimen3*Rb)*EK
1626       if(me.eq.king.or..not.out1file)then
1627        call cartprint
1628        call intout
1629       endif
1630 #ifdef MPI
1631       tt0=MPI_Wtime()
1632 #else
1633       tt0=tcpu()
1634 #endif
1635       call zerograd
1636       call etotal(potEcomp)
1637       if (large) call enerprint(potEcomp)
1638 #ifdef TIMING_ENE
1639 #ifdef MPI
1640       t_etotal=t_etotal+MPI_Wtime()-tt0
1641 #else
1642       t_etotal=t_etotal+tcpu()-tt0
1643 #endif
1644 #endif
1645       potE=potEcomp(0)
1646       call cartgrad
1647       call lagrangian
1648       call max_accel
1649       if (amax*d_time .gt. dvmax) then
1650         d_time=d_time*dvmax/amax
1651         if(me.eq.king.or..not.out1file) write (iout,*)
1652      &   "Time step reduced to",d_time,
1653      &   " because of too large initial acceleration."
1654       endif
1655       if(me.eq.king.or..not.out1file)then
1656        write(iout,*) "Potential energy and its components"
1657        call enerprint(potEcomp)
1658 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1659       endif
1660       potE=potEcomp(0)-potEcomp(20)
1661       totE=EK+potE
1662       itime=0
1663       if (ntwe.ne.0) call statout(itime)
1664       if(me.eq.king.or..not.out1file)
1665      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1666      &   " Kinetic energy",EK," Potential energy",potE,
1667      &   " Total energy",totE," Maximum acceleration ",
1668      &   amax
1669       if (large) then
1670         write (iout,*) "Initial coordinates"
1671         do i=1,nres
1672           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1673      &    (c(j,i+nres),j=1,3)
1674         enddo
1675         write (iout,*) "Initial dC"
1676         do i=0,nres
1677           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1678      &    (dc(j,i+nres),j=1,3)
1679         enddo
1680         write (iout,*) "Initial velocities"
1681         write (iout,"(13x,' backbone ',23x,' side chain')")
1682         do i=0,nres
1683           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1684      &    (d_t(j,i+nres),j=1,3)
1685         enddo
1686         write (iout,*) "Initial accelerations"
1687         do i=0,nres
1688 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1689           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1690      &    (d_a(j,i+nres),j=1,3)
1691         enddo
1692       endif
1693       do i=0,2*nres
1694         do j=1,3
1695           dc_old(j,i)=dc(j,i)
1696           d_t_old(j,i)=d_t(j,i)
1697           d_a_old(j,i)=d_a(j,i)
1698         enddo
1699 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1700       enddo
1701       if (RESPA) then
1702 #ifdef MPI
1703         tt0 =MPI_Wtime()
1704 #else
1705         tt0 = tcpu()
1706 #endif
1707         call zerograd
1708         call etotal_short(energia_short)
1709         if (large) call enerprint(potEcomp)
1710 #ifdef TIMING_ENE
1711 #ifdef MPI
1712         t_eshort=t_eshort+MPI_Wtime()-tt0
1713 #else
1714         t_eshort=t_eshort+tcpu()-tt0
1715 #endif
1716 #endif
1717         call cartgrad
1718         call lagrangian
1719         if(.not.out1file .and. large) then
1720           write (iout,*) "energia_long",energia_long(0),
1721      &     " energia_short",energia_short(0),
1722      &     " total",energia_long(0)+energia_short(0)
1723           write (iout,*) "Initial fast-force accelerations"
1724           do i=0,nres
1725             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1726      &      (d_a(j,i+nres),j=1,3)
1727           enddo
1728         endif
1729 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1730         do i=0,2*nres
1731           do j=1,3
1732             d_a_short(j,i)=d_a(j,i)
1733           enddo
1734         enddo
1735 #ifdef MPI
1736         tt0=MPI_Wtime()
1737 #else
1738         tt0=tcpu()
1739 #endif
1740         call zerograd
1741         call etotal_long(energia_long)
1742         if (large) call enerprint(potEcomp)
1743 #ifdef TIMING_ENE
1744 #ifdef MPI
1745         t_elong=t_elong+MPI_Wtime()-tt0
1746 #else
1747         t_elong=t_elong+tcpu()-tt0
1748 #endif
1749 #endif
1750         call cartgrad
1751         call lagrangian
1752         if(.not.out1file .and. large) then
1753           write (iout,*) "energia_long",energia_long(0)
1754           write (iout,*) "Initial slow-force accelerations"
1755           do i=0,nres
1756             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1757      &      (d_a(j,i+nres),j=1,3)
1758           enddo
1759         endif
1760 #ifdef MPI
1761         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1762 #else
1763         t_enegrad=t_enegrad+tcpu()-tt0
1764 #endif
1765       endif
1766       return
1767       end
1768 c-----------------------------------------------------------
1769       subroutine random_vel
1770       implicit real*8 (a-h,o-z)
1771       include 'DIMENSIONS'
1772       include 'COMMON.CONTROL'
1773       include 'COMMON.VAR'
1774       include 'COMMON.MD'
1775 #ifndef LANG0
1776       include 'COMMON.LANGEVIN'
1777 #else
1778       include 'COMMON.LANGEVIN.lang0'
1779 #endif
1780       include 'COMMON.CHAIN'
1781       include 'COMMON.DERIV'
1782       include 'COMMON.GEO'
1783       include 'COMMON.LOCAL'
1784       include 'COMMON.INTERACT'
1785       include 'COMMON.IOUNITS'
1786       include 'COMMON.NAMES'
1787       include 'COMMON.TIME1'
1788       double precision xv,sigv,lowb,highb
1789 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1790 c First generate velocities in the eigenspace of the G matrix
1791 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1792 c      call flush(iout)
1793       xv=0.0d0
1794       ii=0
1795       do i=1,dimen
1796         do k=1,3
1797           ii=ii+1
1798           sigv=dsqrt((Rb*t_bath)/geigen(i))
1799           lowb=-5*sigv
1800           highb=5*sigv
1801           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1802 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1803 c     &      " d_t_work_new",d_t_work_new(ii)
1804         enddo
1805       enddo
1806 c diagnostics
1807 c      Ek1=0.0d0
1808 c      ii=0
1809 c      do i=1,dimen
1810 c        do k=1,3
1811 c          ii=ii+1
1812 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1813 c        enddo
1814 c      enddo
1815 c      write (iout,*) "Ek from eigenvectors",Ek1
1816 c end diagnostics
1817 c Transform velocities to UNRES coordinate space
1818       do k=0,2
1819         do i=1,dimen
1820           ind=(i-1)*3+k+1
1821           d_t_work(ind)=0.0d0
1822           do j=1,dimen
1823             d_t_work(ind)=d_t_work(ind)
1824      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1825           enddo
1826 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1827 c          call flush(iout)
1828         enddo
1829       enddo
1830 c Transfer to the d_t vector
1831       do j=1,3
1832         d_t(j,0)=d_t_work(j)
1833       enddo
1834       ind=3
1835       do i=nnt,nct-1
1836         do j=1,3
1837           ind=ind+1
1838           d_t(j,i)=d_t_work(ind)
1839         enddo
1840       enddo
1841       do i=nnt,nct
1842         if (itype(i).ne.10 .and. itype(i).ne.21) then
1843           do j=1,3
1844             ind=ind+1
1845             d_t(j,i+nres)=d_t_work(ind)
1846           enddo
1847         endif
1848       enddo
1849 c      call kinetic(EK)
1850 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1851 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1852 c      call flush(iout)
1853       return
1854       end
1855 #ifndef LANG0
1856 c-----------------------------------------------------------
1857       subroutine sd_verlet_p_setup
1858 c Sets up the parameters of stochastic Verlet algorithm
1859       implicit real*8 (a-h,o-z)
1860       include 'DIMENSIONS'
1861 #ifdef MPI
1862       include 'mpif.h'
1863 #endif
1864       include 'COMMON.CONTROL'
1865       include 'COMMON.VAR'
1866       include 'COMMON.MD'
1867 #ifndef LANG0
1868       include 'COMMON.LANGEVIN'
1869 #else
1870       include 'COMMON.LANGEVIN.lang0'
1871 #endif
1872       include 'COMMON.CHAIN'
1873       include 'COMMON.DERIV'
1874       include 'COMMON.GEO'
1875       include 'COMMON.LOCAL'
1876       include 'COMMON.INTERACT'
1877       include 'COMMON.IOUNITS'
1878       include 'COMMON.NAMES'
1879       include 'COMMON.TIME1'
1880       double precision emgdt(MAXRES6),
1881      & pterm,vterm,rho,rhoc,vsig,
1882      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1883      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1884      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1885       logical lprn /.false./
1886       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1887       double precision ktm
1888 #ifdef MPI
1889       tt0 = MPI_Wtime()
1890 #else
1891       tt0 = tcpu()
1892 #endif
1893 c
1894 c AL 8/17/04 Code adapted from tinker
1895 c
1896 c Get the frictional and random terms for stochastic dynamics in the
1897 c eigenspace of mass-scaled UNRES friction matrix
1898 c
1899       do i = 1, dimen
1900             gdt = fricgam(i) * d_time
1901 c
1902 c Stochastic dynamics reduces to simple MD for zero friction
1903 c
1904             if (gdt .le. zero) then
1905                pfric_vec(i) = 1.0d0
1906                vfric_vec(i) = d_time
1907                afric_vec(i) = 0.5d0 * d_time * d_time
1908                prand_vec(i) = 0.0d0
1909                vrand_vec1(i) = 0.0d0
1910                vrand_vec2(i) = 0.0d0
1911 c
1912 c Analytical expressions when friction coefficient is large
1913 c
1914             else
1915                if (gdt .ge. gdt_radius) then
1916                   egdt = dexp(-gdt)
1917                   pfric_vec(i) = egdt
1918                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1919                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1920                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1921                   vterm = 1.0d0 - egdt**2
1922                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1923 c
1924 c Use series expansions when friction coefficient is small
1925 c
1926                else
1927                   gdt2 = gdt * gdt
1928                   gdt3 = gdt * gdt2
1929                   gdt4 = gdt2 * gdt2
1930                   gdt5 = gdt2 * gdt3
1931                   gdt6 = gdt3 * gdt3
1932                   gdt7 = gdt3 * gdt4
1933                   gdt8 = gdt4 * gdt4
1934                   gdt9 = gdt4 * gdt5
1935                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1936      &                          - gdt5/120.0d0 + gdt6/720.0d0
1937      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
1938      &                          - gdt9/362880.0d0) / fricgam(i)**2
1939                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1940                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1941                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1942      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1943      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1944      &                       + 127.0d0*gdt9/90720.0d0
1945                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1946      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1947      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1948      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1949                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1950      &                       - 17.0d0*gdt2/1280.0d0
1951      &                       + 17.0d0*gdt3/6144.0d0
1952      &                       + 40967.0d0*gdt4/34406400.0d0
1953      &                       - 57203.0d0*gdt5/275251200.0d0
1954      &                       - 1429487.0d0*gdt6/13212057600.0d0)
1955                end if
1956 c
1957 c Compute the scaling factors of random terms for the nonzero friction case
1958 c
1959                ktm = 0.5d0*d_time/fricgam(i)
1960                psig = dsqrt(ktm*pterm) / fricgam(i)
1961                vsig = dsqrt(ktm*vterm)
1962                rhoc = dsqrt(1.0d0 - rho*rho)
1963                prand_vec(i) = psig
1964                vrand_vec1(i) = vsig * rho
1965                vrand_vec2(i) = vsig * rhoc
1966             end if
1967       end do
1968       if (lprn) then
1969       write (iout,*)
1970      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1971      &  " vrand_vec2"
1972       do i=1,dimen
1973         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1974      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1975       enddo
1976       endif
1977 c
1978 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1979 c
1980 #ifndef   LANG0
1981       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1982       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1983       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1984       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1985       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1986       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1987 #endif
1988 #ifdef MPI
1989       t_sdsetup=t_sdsetup+MPI_Wtime()
1990 #else
1991       t_sdsetup=t_sdsetup+tcpu()-tt0
1992 #endif
1993       return
1994       end
1995 c-------------------------------------------------------------
1996       subroutine eigtransf1(n,ndim,ab,d,c)
1997       implicit none
1998       integer n,ndim
1999       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2000       integer i,j,k
2001       do i=1,n
2002         do j=1,n
2003           c(i,j)=0.0d0
2004           do k=1,n
2005             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2006           enddo
2007         enddo
2008       enddo
2009       return
2010       end
2011 c-------------------------------------------------------------
2012       subroutine eigtransf(n,ndim,a,b,d,c)
2013       implicit none
2014       integer n,ndim
2015       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2016       integer i,j,k
2017       do i=1,n
2018         do j=1,n
2019           c(i,j)=0.0d0
2020           do k=1,n
2021             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2022           enddo
2023         enddo
2024       enddo
2025       return
2026       end
2027 c-------------------------------------------------------------
2028       subroutine sd_verlet1
2029 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2030       implicit real*8 (a-h,o-z)
2031       include 'DIMENSIONS'
2032       include 'COMMON.CONTROL'
2033       include 'COMMON.VAR'
2034       include 'COMMON.MD'
2035 #ifndef LANG0
2036       include 'COMMON.LANGEVIN'
2037 #else
2038       include 'COMMON.LANGEVIN.lang0'
2039 #endif
2040       include 'COMMON.CHAIN'
2041       include 'COMMON.DERIV'
2042       include 'COMMON.GEO'
2043       include 'COMMON.LOCAL'
2044       include 'COMMON.INTERACT'
2045       include 'COMMON.IOUNITS'
2046       include 'COMMON.NAMES'
2047       double precision stochforcvec(MAXRES6)
2048       common /stochcalc/ stochforcvec
2049       logical lprn /.false./
2050
2051 c      write (iout,*) "dc_old"
2052 c      do i=0,nres
2053 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2054 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2055 c      enddo
2056       do j=1,3
2057         dc_work(j)=dc_old(j,0)
2058         d_t_work(j)=d_t_old(j,0)
2059         d_a_work(j)=d_a_old(j,0)
2060       enddo
2061       ind=3
2062       do i=nnt,nct-1
2063         do j=1,3
2064           dc_work(ind+j)=dc_old(j,i)
2065           d_t_work(ind+j)=d_t_old(j,i)
2066           d_a_work(ind+j)=d_a_old(j,i)
2067         enddo
2068         ind=ind+3
2069       enddo
2070       do i=nnt,nct
2071         if (itype(i).ne.10 .and. itype(i).ne.21) then
2072           do j=1,3
2073             dc_work(ind+j)=dc_old(j,i+nres)
2074             d_t_work(ind+j)=d_t_old(j,i+nres)
2075             d_a_work(ind+j)=d_a_old(j,i+nres)
2076           enddo
2077           ind=ind+3
2078         endif
2079       enddo
2080 #ifndef LANG0
2081       if (lprn) then
2082       write (iout,*)
2083      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2084      &  " vrand_mat2"
2085       do i=1,dimen
2086         do j=1,dimen
2087           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2088      &      vfric_mat(i,j),afric_mat(i,j),
2089      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2090         enddo
2091       enddo
2092       endif
2093       do i=1,dimen
2094         ddt1=0.0d0
2095         ddt2=0.0d0
2096         do j=1,dimen
2097           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2098      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2099           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2100           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2101         enddo
2102         d_t_work_new(i)=ddt1+0.5d0*ddt2
2103         d_t_work(i)=ddt1+ddt2
2104       enddo
2105 #endif
2106       do j=1,3
2107         dc(j,0)=dc_work(j)
2108         d_t(j,0)=d_t_work(j)
2109       enddo
2110       ind=3
2111       do i=nnt,nct-1
2112         do j=1,3
2113           dc(j,i)=dc_work(ind+j)
2114           d_t(j,i)=d_t_work(ind+j)
2115         enddo
2116         ind=ind+3
2117       enddo
2118       do i=nnt,nct
2119         if (itype(i).ne.10) then
2120           inres=i+nres
2121           do j=1,3
2122             dc(j,inres)=dc_work(ind+j)
2123             d_t(j,inres)=d_t_work(ind+j)
2124           enddo
2125           ind=ind+3
2126         endif
2127       enddo
2128       return
2129       end
2130 c--------------------------------------------------------------------------
2131       subroutine sd_verlet2
2132 c  Calculating the adjusted velocities for accelerations
2133       implicit real*8 (a-h,o-z)
2134       include 'DIMENSIONS'
2135       include 'COMMON.CONTROL'
2136       include 'COMMON.VAR'
2137       include 'COMMON.MD'
2138 #ifndef LANG0
2139       include 'COMMON.LANGEVIN'
2140 #else
2141       include 'COMMON.LANGEVIN.lang0'
2142 #endif
2143       include 'COMMON.CHAIN'
2144       include 'COMMON.DERIV'
2145       include 'COMMON.GEO'
2146       include 'COMMON.LOCAL'
2147       include 'COMMON.INTERACT'
2148       include 'COMMON.IOUNITS'
2149       include 'COMMON.NAMES'
2150       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2151       common /stochcalc/ stochforcvec
2152 c
2153 c Compute the stochastic forces which contribute to velocity change
2154 c
2155       call stochastic_force(stochforcvecV)
2156
2157 #ifndef LANG0
2158       do i=1,dimen
2159         ddt1=0.0d0
2160         ddt2=0.0d0
2161         do j=1,dimen
2162           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2163           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2164      &     vrand_mat2(i,j)*stochforcvecV(j)
2165         enddo
2166         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2167       enddo
2168 #endif
2169       do j=1,3
2170         d_t(j,0)=d_t_work(j)
2171       enddo
2172       ind=3
2173       do i=nnt,nct-1
2174         do j=1,3
2175           d_t(j,i)=d_t_work(ind+j)
2176         enddo
2177         ind=ind+3
2178       enddo
2179       do i=nnt,nct
2180         if (itype(i).ne.10 .and. itype(i).ne.21) then
2181           inres=i+nres
2182           do j=1,3
2183             d_t(j,inres)=d_t_work(ind+j)
2184           enddo
2185           ind=ind+3
2186         endif
2187       enddo
2188       return
2189       end
2190 c-----------------------------------------------------------
2191       subroutine sd_verlet_ciccotti_setup
2192 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2193 c version
2194       implicit real*8 (a-h,o-z)
2195       include 'DIMENSIONS'
2196 #ifdef MPI
2197       include 'mpif.h'
2198 #endif
2199       include 'COMMON.CONTROL'
2200       include 'COMMON.VAR'
2201       include 'COMMON.MD'
2202 #ifndef LANG0
2203       include 'COMMON.LANGEVIN'
2204 #else
2205       include 'COMMON.LANGEVIN.lang0'
2206 #endif
2207       include 'COMMON.CHAIN'
2208       include 'COMMON.DERIV'
2209       include 'COMMON.GEO'
2210       include 'COMMON.LOCAL'
2211       include 'COMMON.INTERACT'
2212       include 'COMMON.IOUNITS'
2213       include 'COMMON.NAMES'
2214       include 'COMMON.TIME1'
2215       double precision emgdt(MAXRES6),
2216      & pterm,vterm,rho,rhoc,vsig,
2217      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2218      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2219      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2220       logical lprn /.false./
2221       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2222       double precision ktm
2223 #ifdef MPI
2224       tt0 = MPI_Wtime()
2225 #else
2226       tt0 = tcpu()
2227 #endif
2228 c
2229 c AL 8/17/04 Code adapted from tinker
2230 c
2231 c Get the frictional and random terms for stochastic dynamics in the
2232 c eigenspace of mass-scaled UNRES friction matrix
2233 c
2234       do i = 1, dimen
2235             write (iout,*) "i",i," fricgam",fricgam(i)
2236             gdt = fricgam(i) * d_time
2237 c
2238 c Stochastic dynamics reduces to simple MD for zero friction
2239 c
2240             if (gdt .le. zero) then
2241                pfric_vec(i) = 1.0d0
2242                vfric_vec(i) = d_time
2243                afric_vec(i) = 0.5d0*d_time*d_time
2244                prand_vec(i) = afric_vec(i)
2245                vrand_vec2(i) = vfric_vec(i)
2246 c
2247 c Analytical expressions when friction coefficient is large
2248 c
2249             else
2250                egdt = dexp(-gdt)
2251                pfric_vec(i) = egdt
2252                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2253                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2254                prand_vec(i) = afric_vec(i)
2255                vrand_vec2(i) = vfric_vec(i)
2256 c
2257 c Compute the scaling factors of random terms for the nonzero friction case
2258 c
2259 c               ktm = 0.5d0*d_time/fricgam(i)
2260 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2261 c               vsig = dsqrt(ktm*vterm)
2262 c               prand_vec(i) = psig*afric_vec(i)
2263 c               vrand_vec2(i) = vsig*vfric_vec(i)
2264             end if
2265       end do
2266       if (lprn) then
2267       write (iout,*)
2268      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2269      &  " vrand_vec2"
2270       do i=1,dimen
2271         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2272      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2273       enddo
2274       endif
2275 c
2276 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2277 c
2278       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2279       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2280       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2281       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2282       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2283 #ifdef MPI
2284       t_sdsetup=t_sdsetup+MPI_Wtime()
2285 #else
2286       t_sdsetup=t_sdsetup+tcpu()-tt0
2287 #endif
2288       return
2289       end
2290 c-------------------------------------------------------------
2291       subroutine sd_verlet1_ciccotti
2292 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2293       implicit real*8 (a-h,o-z)
2294       include 'DIMENSIONS'
2295 #ifdef MPI
2296       include 'mpif.h'
2297 #endif
2298       include 'COMMON.CONTROL'
2299       include 'COMMON.VAR'
2300       include 'COMMON.MD'
2301 #ifndef LANG0
2302       include 'COMMON.LANGEVIN'
2303 #else
2304       include 'COMMON.LANGEVIN.lang0'
2305 #endif
2306       include 'COMMON.CHAIN'
2307       include 'COMMON.DERIV'
2308       include 'COMMON.GEO'
2309       include 'COMMON.LOCAL'
2310       include 'COMMON.INTERACT'
2311       include 'COMMON.IOUNITS'
2312       include 'COMMON.NAMES'
2313       double precision stochforcvec(MAXRES6)
2314       common /stochcalc/ stochforcvec
2315       logical lprn /.false./
2316
2317 c      write (iout,*) "dc_old"
2318 c      do i=0,nres
2319 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2320 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2321 c      enddo
2322       do j=1,3
2323         dc_work(j)=dc_old(j,0)
2324         d_t_work(j)=d_t_old(j,0)
2325         d_a_work(j)=d_a_old(j,0)
2326       enddo
2327       ind=3
2328       do i=nnt,nct-1
2329         do j=1,3
2330           dc_work(ind+j)=dc_old(j,i)
2331           d_t_work(ind+j)=d_t_old(j,i)
2332           d_a_work(ind+j)=d_a_old(j,i)
2333         enddo
2334         ind=ind+3
2335       enddo
2336       do i=nnt,nct
2337         if (itype(i).ne.10 .and. itype(i).ne.21) then
2338           do j=1,3
2339             dc_work(ind+j)=dc_old(j,i+nres)
2340             d_t_work(ind+j)=d_t_old(j,i+nres)
2341             d_a_work(ind+j)=d_a_old(j,i+nres)
2342           enddo
2343           ind=ind+3
2344         endif
2345       enddo
2346
2347 #ifndef LANG0
2348       if (lprn) then
2349       write (iout,*)
2350      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2351      &  " vrand_mat2"
2352       do i=1,dimen
2353         do j=1,dimen
2354           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2355      &      vfric_mat(i,j),afric_mat(i,j),
2356      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2357         enddo
2358       enddo
2359       endif
2360       do i=1,dimen
2361         ddt1=0.0d0
2362         ddt2=0.0d0
2363         do j=1,dimen
2364           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2365      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2366           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2367           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2368         enddo
2369         d_t_work_new(i)=ddt1+0.5d0*ddt2
2370         d_t_work(i)=ddt1+ddt2
2371       enddo
2372 #endif
2373       do j=1,3
2374         dc(j,0)=dc_work(j)
2375         d_t(j,0)=d_t_work(j)
2376       enddo
2377       ind=3
2378       do i=nnt,nct-1
2379         do j=1,3
2380           dc(j,i)=dc_work(ind+j)
2381           d_t(j,i)=d_t_work(ind+j)
2382         enddo
2383         ind=ind+3
2384       enddo
2385       do i=nnt,nct
2386         if (itype(i).ne.10 .and. itype(i).ne.21) then
2387           inres=i+nres
2388           do j=1,3
2389             dc(j,inres)=dc_work(ind+j)
2390             d_t(j,inres)=d_t_work(ind+j)
2391           enddo
2392           ind=ind+3
2393         endif
2394       enddo
2395       return
2396       end
2397 c--------------------------------------------------------------------------
2398       subroutine sd_verlet2_ciccotti
2399 c  Calculating the adjusted velocities for accelerations
2400       implicit real*8 (a-h,o-z)
2401       include 'DIMENSIONS'
2402       include 'COMMON.CONTROL'
2403       include 'COMMON.VAR'
2404       include 'COMMON.MD'
2405 #ifndef LANG0
2406       include 'COMMON.LANGEVIN'
2407 #else
2408       include 'COMMON.LANGEVIN.lang0'
2409 #endif
2410       include 'COMMON.CHAIN'
2411       include 'COMMON.DERIV'
2412       include 'COMMON.GEO'
2413       include 'COMMON.LOCAL'
2414       include 'COMMON.INTERACT'
2415       include 'COMMON.IOUNITS'
2416       include 'COMMON.NAMES'
2417       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2418       common /stochcalc/ stochforcvec
2419 c
2420 c Compute the stochastic forces which contribute to velocity change
2421 c
2422       call stochastic_force(stochforcvecV)
2423 #ifndef LANG0
2424       do i=1,dimen
2425         ddt1=0.0d0
2426         ddt2=0.0d0
2427         do j=1,dimen
2428
2429           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2430 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2431           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2432         enddo
2433         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2434       enddo
2435 #endif
2436       do j=1,3
2437         d_t(j,0)=d_t_work(j)
2438       enddo
2439       ind=3
2440       do i=nnt,nct-1
2441         do j=1,3
2442           d_t(j,i)=d_t_work(ind+j)
2443         enddo
2444         ind=ind+3
2445       enddo
2446       do i=nnt,nct
2447         if (itype(i).ne.10 .and. itype(i).ne.21) then
2448           inres=i+nres
2449           do j=1,3
2450             d_t(j,inres)=d_t_work(ind+j)
2451           enddo
2452           ind=ind+3
2453         endif
2454       enddo
2455       return
2456       end
2457 #endif