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

   1       subroutine MD
   2 c------------------------------------------------
   3 c  The driver for molecular dynamics subroutines
   4 c------------------------------------------------
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7 #ifdef MPI
   8       include "mpif.h"
   9       integer IERROR,ERRCODE
  10 #endif
  11       include 'COMMON.SETUP'
  12       include 'COMMON.CONTROL'
  13       include 'COMMON.VAR'
  14       include 'COMMON.MD'
  15 #ifndef LANG0
  16       include 'COMMON.LANGEVIN'
  17 #else
  18       include 'COMMON.LANGEVIN.lang0'
  19 #endif
  20       include 'COMMON.CHAIN'
  21       include 'COMMON.DERIV'
  22       include 'COMMON.GEO'
  23       include 'COMMON.LOCAL'
  24       include 'COMMON.INTERACT'
  25       include 'COMMON.IOUNITS'
  26       include 'COMMON.NAMES'
  27       include 'COMMON.TIME1'
  28       include 'COMMON.HAIRPIN'
  29       double precision cm(3),L(3),vcm(3)
  30 #ifdef VOUT
  31       double precision v_work(maxres6),v_transf(maxres6)
  32 #endif
  33       integer ilen,rstcount
  34       external ilen
  35       character*50 tytul
  36       common /gucio/ cm
  37       integer itime
  38       logical ovrtim
  39 c
  40 #ifdef MPI
  41       if (ilen(tmpdir).gt.0)
  42      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
  43      &        //liczba(:ilen(liczba))//'.rst')
  44 #else
  45       if (ilen(tmpdir).gt.0)
  46      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
  47 #endif
  48       t_MDsetup=0.0d0
  49       t_langsetup=0.0d0
  50       t_MD=0.0d0
  51       t_enegrad=0.0d0
  52       t_sdsetup=0.0d0
  53       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  54 #ifdef MPI
  55       tt0=MPI_Wtime()
  56 #else
  57       tt0 = tcpu()
  58 #endif
  59 c Determine the inverse of the inertia matrix.
  60       call setup_MD_matrices
  61 c Initialize MD
  62       call init_MD
  63 #ifdef MPI
  64       t_MDsetup = MPI_Wtime()-tt0
  65 #else
  66       t_MDsetup = tcpu()-tt0
  67 #endif
  68       rstcount=0
  69 c   Entering the MD loop
  70 #ifdef MPI
  71       tt0 = MPI_Wtime()
  72 #else
  73       tt0 = tcpu()
  74 #endif
  75       if (lang.eq.2 .or. lang.eq.3) then
  76 #ifndef   LANG0
  77         call setup_fricmat
  78         if (lang.eq.2) then
  79           call sd_verlet_p_setup
  80         else
  81           call sd_verlet_ciccotti_setup
  82         endif
  83         do i=1,dimen
  84           do j=1,dimen
  85             pfric0_mat(i,j,0)=pfric_mat(i,j)
  86             afric0_mat(i,j,0)=afric_mat(i,j)
  87             vfric0_mat(i,j,0)=vfric_mat(i,j)
  88             prand0_mat(i,j,0)=prand_mat(i,j)
  89             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
  90             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
  91           enddo
  92         enddo
  93         flag_stoch(0)=.true.
  94         do i=1,maxflag_stoch
  95           flag_stoch(i)=.false.
  96         enddo
  97 #else
  98         write (iout,*)
  99      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 100 #ifdef MPI
 101         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 102 #endif
 103         stop
 104 #endif
 105       else if (lang.eq.1 .or. lang.eq.4) then
 106         call setup_fricmat
 107       endif
 108 #ifdef MPI
 109       t_langsetup=MPI_Wtime()-tt0
 110       tt0=MPI_Wtime()
 111 #else
 112       t_langsetup=tcpu()-tt0
 113       tt0=tcpu()
 114 #endif
 115       do itime=1,n_timestep
 116         if (ovrtim()) exit
 117         if (large.and. mod(itime,ntwe).eq.0)
 118      &    write (iout,*) "itime",itime
 119         rstcount=rstcount+1
 120         if (lang.gt.0 .and. surfarea .and.
 121      &      mod(itime,reset_fricmat).eq.0) then
 122           if (lang.eq.2 .or. lang.eq.3) then
 123 #ifndef LANG0
 124             call setup_fricmat
 125             if (lang.eq.2) then
 126               call sd_verlet_p_setup
 127             else
 128               call sd_verlet_ciccotti_setup
 129             endif
 130             do i=1,dimen
 131               do j=1,dimen
 132                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 133                 afric0_mat(i,j,0)=afric_mat(i,j)
 134                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 135                 prand0_mat(i,j,0)=prand_mat(i,j)
 136                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 137                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 138               enddo
 139             enddo
 140             flag_stoch(0)=.true.
 141             do i=1,maxflag_stoch
 142               flag_stoch(i)=.false.
 143             enddo
 144 #endif
 145           else if (lang.eq.1 .or. lang.eq.4) then
 146             call setup_fricmat
 147           endif
 148           write (iout,'(a,i10)')
 149      &      "Friction matrix reset based on surface area, itime",itime
 150         endif
 151         if (reset_vel .and. tbf .and. lang.eq.0
 152      &      .and. mod(itime,count_reset_vel).eq.0) then
 153           call random_vel
 154           write(iout,'(a,f20.2)')
 155      &     "Velocities reset to random values, time",totT
 156           do i=0,2*nres
 157             do j=1,3
 158               d_t_old(j,i)=d_t(j,i)
 159             enddo
 160           enddo
 161         endif
 162         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 163           call inertia_tensor
 164           call vcm_vel(vcm)
 165           do j=1,3
 166              d_t(j,0)=d_t(j,0)-vcm(j)
 167           enddo
 168           call kinetic(EK)
 169           kinetic_T=2.0d0/(dimen3*Rb)*EK
 170           scalfac=dsqrt(T_bath/kinetic_T)
 171           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 172           do i=0,2*nres
 173             do j=1,3
 174               d_t_old(j,i)=scalfac*d_t(j,i)
 175             enddo
 176           enddo
 177         endif
 178         if (lang.ne.4) then
 179           if (RESPA) then
 180 c Time-reversible RESPA algorithm
 181 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 182             call RESPA_step(itime)
 183           else
 184 c Variable time step algorithm.
 185             call velverlet_step(itime)
 186           endif
 187         else
 188 #ifdef BROWN
 189           call brown_step(itime)
 190 #else
 191           print *,"Brown dynamics not here!"
 192 #ifdef MPI
 193           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 194 #endif
 195           stop
 196 #endif
 197         endif
 198         if (ntwe.ne.0) then
 199          if (mod(itime,ntwe).eq.0) call statout(itime)
 200 #ifdef VOUT
 201         do j=1,3
 202           v_work(j)=d_t(j,0)
 203         enddo
 204         ind=3
 205         do i=nnt,nct-1
 206           do j=1,3
 207             ind=ind+1
 208             v_work(ind)=d_t(j,i)
 209           enddo
 210         enddo
 211         do i=nnt,nct
 212           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 213             do j=1,3
 214               ind=ind+1
 215               v_work(ind)=d_t(j,i+nres)
 216             enddo
 217           endif
 218         enddo
 219
 220         write (66,'(80f10.5)')
 221      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 222         do i=1,ind
 223           v_transf(i)=0.0d0
 224           do j=1,ind
 225             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 226           enddo
 227            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 228         enddo
 229         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 230 #endif
 231         endif
 232         if (mod(itime,ntwx).eq.0) then
 233           write (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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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.ntyp1) 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 c              write (iout,*) i,dacc
1264               if (dacc.gt.amax) amax=dacc
1265             endif
1266           enddo
1267         endif
1268       enddo
1269 c Side chains
1270       do j=1,3
1271 c        accel(j)=aux(j)
1272         accel_old(j)=d_a_old(j,0)
1273         accel(j)=d_a(j,0)
1274       enddo
1275       if (nnt.eq.2) then
1276         do j=1,3
1277           accel_old(j)=accel_old(j)+d_a_old(j,1)
1278           accel(j)=accel(j)+d_a(j,1)
1279         enddo
1280       endif
1281       do i=nnt,nct
1282         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1283           do j=1,3
1284 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1285             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1286             accel(j)=accel(j)+d_a(j,i+nres)
1287           enddo
1288         endif
1289         do j=1,3
1290 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1291           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1292             dacc=dabs(accel(j)-accel_old(j))
1293 c            write (iout,*) "side-chain",i,dacc
1294             if (dacc.gt.amax) amax=dacc
1295           endif
1296         enddo
1297         do j=1,3
1298           accel_old(j)=accel_old(j)+d_a_old(j,i)
1299           accel(j)=accel(j)+d_a(j,i)
1300 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1301         enddo
1302       enddo
1303       return
1304       end
1305 c---------------------------------------------------------------------
1306       subroutine predict_edrift(epdrift)
1307 c
1308 c Predict the drift of the potential energy
1309 c
1310       implicit real*8 (a-h,o-z)
1311       include 'DIMENSIONS'
1312       include 'COMMON.CONTROL'
1313       include 'COMMON.VAR'
1314       include 'COMMON.MD'
1315       include 'COMMON.CHAIN'
1316       include 'COMMON.DERIV'
1317       include 'COMMON.GEO'
1318       include 'COMMON.LOCAL'
1319       include 'COMMON.INTERACT'
1320       include 'COMMON.IOUNITS'
1321       include 'COMMON.MUCA'
1322       double precision epdrift,epdriftij
1323 c Drift of the potential energy
1324       epdrift=0.0d0
1325       do i=nnt,nct
1326 c Backbone
1327         if (i.lt.nct) then
1328           do j=1,3
1329             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1330             if (lmuca) epdriftij=epdriftij*factor
1331 c            write (iout,*) "back",i,j,epdriftij
1332             if (epdriftij.gt.epdrift) epdrift=epdriftij
1333           enddo
1334         endif
1335 c Side chains
1336         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1337           do j=1,3
1338             epdriftij=
1339      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1340             if (lmuca) epdriftij=epdriftij*factor
1341 c            write (iout,*) "side",i,j,epdriftij
1342             if (epdriftij.gt.epdrift) epdrift=epdriftij
1343           enddo
1344         endif
1345       enddo
1346       epdrift=0.5d0*epdrift*d_time*d_time
1347 c      write (iout,*) "epdrift",epdrift
1348       return
1349       end
1350 c-----------------------------------------------------------------------
1351       subroutine verlet_bath
1352 c
1353 c  Coupling to the thermostat by using the Berendsen algorithm
1354 c
1355       implicit real*8 (a-h,o-z)
1356       include 'DIMENSIONS'
1357       include 'COMMON.CONTROL'
1358       include 'COMMON.VAR'
1359       include 'COMMON.MD'
1360       include 'COMMON.CHAIN'
1361       include 'COMMON.DERIV'
1362       include 'COMMON.GEO'
1363       include 'COMMON.LOCAL'
1364       include 'COMMON.INTERACT'
1365       include 'COMMON.IOUNITS'
1366       include 'COMMON.NAMES'
1367       double precision T_half,fact
1368 c
1369       T_half=2.0d0/(dimen3*Rb)*EK
1370       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1371 c      write(iout,*) "T_half", T_half
1372 c      write(iout,*) "EK", EK
1373 c      write(iout,*) "fact", fact
1374       do j=1,3
1375         d_t(j,0)=fact*d_t(j,0)
1376       enddo
1377       do i=nnt,nct-1
1378         do j=1,3
1379           d_t(j,i)=fact*d_t(j,i)
1380         enddo
1381       enddo
1382       do i=nnt,nct
1383         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1384           inres=i+nres
1385           do j=1,3
1386             d_t(j,inres)=fact*d_t(j,inres)
1387           enddo
1388         endif
1389       enddo
1390       return
1391       end
1392 c---------------------------------------------------------
1393       subroutine init_MD
1394 c  Set up the initial conditions of a MD simulation
1395       implicit real*8 (a-h,o-z)
1396       include 'DIMENSIONS'
1397 #ifdef MP
1398       include 'mpif.h'
1399       character*16 form
1400       integer IERROR,ERRCODE
1401 #endif
1402       include 'COMMON.SETUP'
1403       include 'COMMON.CONTROL'
1404       include 'COMMON.VAR'
1405       include 'COMMON.MD'
1406 #ifndef LANG0
1407       include 'COMMON.LANGEVIN'
1408 #else
1409       include 'COMMON.LANGEVIN.lang0'
1410 #endif
1411       include 'COMMON.CHAIN'
1412       include 'COMMON.DERIV'
1413       include 'COMMON.GEO'
1414       include 'COMMON.LOCAL'
1415       include 'COMMON.INTERACT'
1416       include 'COMMON.IOUNITS'
1417       include 'COMMON.NAMES'
1418       include 'COMMON.REMD'
1419       real*8 energia_long(0:n_ene),
1420      &  energia_short(0:n_ene),vcm(3),incr(3)
1421       double precision cm(3),L(3),xv,sigv,lowb,highb
1422       double precision varia(maxvar)
1423       character*256 qstr
1424       integer ilen
1425       external ilen
1426       character*50 tytul
1427       logical file_exist
1428       common /gucio/ cm
1429       d_time0=d_time
1430 c      write(iout,*) "d_time", d_time
1431 c Compute the standard deviations of stochastic forces for Langevin dynamics
1432 c if the friction coefficients do not depend on surface area
1433       if (lang.gt.0 .and. .not.surfarea) then
1434         do i=nnt,nct-1
1435           stdforcp(i)=stdfp*dsqrt(gamp)
1436         enddo
1437         do i=nnt,nct
1438           stdforcsc(i)=stdfsc(iabs(itype(i)))
1439      &                *dsqrt(gamsc(iabs(itype(i))))
1440         enddo
1441       endif
1442 c Open the pdb file for snapshotshots
1443 #ifdef MPI
1444       if(mdpdb) then
1445         if (ilen(tmpdir).gt.0)
1446      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1447      &      liczba(:ilen(liczba))//".pdb")
1448         open(ipdb,
1449      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1450      &  //".pdb")
1451       else
1452 #ifdef NOXDR
1453         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1454      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1455      &      liczba(:ilen(liczba))//".x")
1456         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1457      &  //".x"
1458 #else
1459         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1460      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1461      &      liczba(:ilen(liczba))//".cx")
1462         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1463      &  //".cx"
1464 #endif
1465       endif
1466 #else
1467       if(mdpdb) then
1468          if (ilen(tmpdir).gt.0)
1469      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1470          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1471       else
1472          if (ilen(tmpdir).gt.0)
1473      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1474          cartname=prefix(:ilen(prefix))//"_MD.cx"
1475       endif
1476 #endif
1477       if (usampl) then
1478         write (qstr,'(256(1h ))')
1479         ipos=1
1480         do i=1,nfrag
1481           iq = qinfrag(i,iset)*10
1482           iw = wfrag(i,iset)/100
1483           if (iw.gt.0) then
1484             if(me.eq.king.or..not.out1file)
1485      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1486             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1487             ipos=ipos+7
1488           endif
1489         enddo
1490         do i=1,npair
1491           iq = qinpair(i,iset)*10
1492           iw = wpair(i,iset)/100
1493           if (iw.gt.0) then
1494             if(me.eq.king.or..not.out1file)
1495      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1496             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1497             ipos=ipos+7
1498           endif
1499         enddo
1500 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1501 #ifdef NOXDR
1502 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1503 #else
1504 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1505 #endif
1506 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1507       endif
1508       icg=1
1509       if (rest) then
1510        if (restart1file) then
1511          if (me.eq.king)
1512      &     inquire(file=mremd_rst_name,exist=file_exist)
1513            write (*,*) me," Before broadcast: file_exist",file_exist
1514          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1515      &          IERR)
1516          write (*,*) me," After broadcast: file_exist",file_exist
1517 c        inquire(file=mremd_rst_name,exist=file_exist)
1518         if(me.eq.king.or..not.out1file)
1519      &   write(iout,*) "Initial state read by master and distributed"
1520        else
1521          if (ilen(tmpdir).gt.0)
1522      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1523      &      //liczba(:ilen(liczba))//'.rst')
1524         inquire(file=rest2name,exist=file_exist)
1525        endif
1526        if(file_exist) then
1527          if(.not.restart1file) then
1528            if(me.eq.king.or..not.out1file)
1529      &      write(iout,*) "Initial state will be read from file ",
1530      &      rest2name(:ilen(rest2name))
1531            call readrst
1532          endif
1533          call rescale_weights(t_bath)
1534        else
1535         if(me.eq.king.or..not.out1file)then
1536          if (restart1file) then
1537           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1538      &       " does not exist"
1539          else
1540           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1541      &       " does not exist"
1542          endif
1543          write(iout,*) "Initial velocities randomly generated"
1544         endif
1545         call random_vel
1546         totT=0.0d0
1547        endif
1548       else
1549 c Generate initial velocities
1550         if(me.eq.king.or..not.out1file)
1551      &   write(iout,*) "Initial velocities randomly generated"
1552         call random_vel
1553         totT=0.0d0
1554       endif
1555 c      rest2name = prefix(:ilen(prefix))//'.rst'
1556       if(me.eq.king.or..not.out1file)then
1557        write (iout,*) "Initial velocities"
1558        do i=0,nres
1559          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1560      &   (d_t(j,i+nres),j=1,3)
1561        enddo
1562 c  Zeroing the total angular momentum of the system
1563        write(iout,*) "Calling the zero-angular
1564      & momentum subroutine"
1565       endif
1566       call inertia_tensor
1567 c  Getting the potential energy and forces and velocities and accelerations
1568       call vcm_vel(vcm)
1569 c      write (iout,*) "velocity of the center of the mass:"
1570 c      write (iout,*) (vcm(j),j=1,3)
1571       do j=1,3
1572         d_t(j,0)=d_t(j,0)-vcm(j)
1573       enddo
1574 c Removing the velocity of the center of mass
1575       call vcm_vel(vcm)
1576       if(me.eq.king.or..not.out1file)then
1577        write (iout,*) "vcm right after adjustment:"
1578        write (iout,*) (vcm(j),j=1,3)
1579       endif
1580       if (.not.rest) then
1581          call chainbuild
1582          if(iranconf.ne.0) then
1583           if (overlapsc) then
1584            print *, 'Calling OVERLAP_SC'
1585            call overlap_sc(fail)
1586           endif
1587
1588           if (searchsc) then
1589            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1590            print *,'SC_move',nft_sc,etot
1591            if(me.eq.king.or..not.out1file)
1592      &      write(iout,*) 'SC_move',nft_sc,etot
1593           endif
1594
1595           if(dccart)then
1596            print *, 'Calling MINIM_DC'
1597            call minim_dc(etot,iretcode,nfun)
1598           else
1599            call geom_to_var(nvar,varia)
1600            print *,'Calling MINIMIZE.'
1601            call minimize(etot,varia,iretcode,nfun)
1602            call var_to_geom(nvar,varia)
1603           endif
1604           if(me.eq.king.or..not.out1file)
1605      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1606          endif
1607       endif
1608       call chainbuild_cart
1609       call kinetic(EK)
1610       if (tbf) then
1611         call verlet_bath
1612       endif
1613       kinetic_T=2.0d0/(dimen3*Rb)*EK
1614       if(me.eq.king.or..not.out1file)then
1615        call cartprint
1616        call intout
1617       endif
1618 #ifdef MPI
1619       tt0=MPI_Wtime()
1620 #else
1621       tt0=tcpu()
1622 #endif
1623       call zerograd
1624       call etotal(potEcomp)
1625       if (large) call enerprint(potEcomp)
1626 #ifdef TIMING_ENE
1627 #ifdef MPI
1628       t_etotal=t_etotal+MPI_Wtime()-tt0
1629 #else
1630       t_etotal=t_etotal+tcpu()-tt0
1631 #endif
1632 #endif
1633       potE=potEcomp(0)
1634       call cartgrad
1635       call lagrangian
1636       call max_accel
1637       if (amax*d_time .gt. dvmax) then
1638         d_time=d_time*dvmax/amax
1639         if(me.eq.king.or..not.out1file) write (iout,*)
1640      &   "Time step reduced to",d_time,
1641      &   " because of too large initial acceleration."
1642       endif
1643       if(me.eq.king.or..not.out1file)then
1644        write(iout,*) "Potential energy and its components"
1645        call enerprint(potEcomp)
1646 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1647       endif
1648       potE=potEcomp(0)-potEcomp(20)
1649       totE=EK+potE
1650       itime=0
1651       if (ntwe.ne.0) call statout(itime)
1652       if(me.eq.king.or..not.out1file)
1653      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1654      &   " Kinetic energy",EK," Potential energy",potE,
1655      &   " Total energy",totE," Maximum acceleration ",
1656      &   amax
1657       if (large) then
1658         write (iout,*) "Initial coordinates"
1659         do i=1,nres
1660           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1661      &    (c(j,i+nres),j=1,3)
1662         enddo
1663         write (iout,*) "Initial dC"
1664         do i=0,nres
1665           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1666      &    (dc(j,i+nres),j=1,3)
1667         enddo
1668         write (iout,*) "Initial velocities"
1669         write (iout,"(13x,' backbone ',23x,' side chain')")
1670         do i=0,nres
1671           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1672      &    (d_t(j,i+nres),j=1,3)
1673         enddo
1674         write (iout,*) "Initial accelerations"
1675         do i=0,nres
1676 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1677           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1678      &    (d_a(j,i+nres),j=1,3)
1679         enddo
1680       endif
1681       do i=0,2*nres
1682         do j=1,3
1683           dc_old(j,i)=dc(j,i)
1684           d_t_old(j,i)=d_t(j,i)
1685           d_a_old(j,i)=d_a(j,i)
1686         enddo
1687 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1688       enddo
1689       if (RESPA) then
1690 #ifdef MPI
1691         tt0 =MPI_Wtime()
1692 #else
1693         tt0 = tcpu()
1694 #endif
1695         call zerograd
1696         call etotal_short(energia_short)
1697         if (large) call enerprint(potEcomp)
1698 #ifdef TIMING_ENE
1699 #ifdef MPI
1700         t_eshort=t_eshort+MPI_Wtime()-tt0
1701 #else
1702         t_eshort=t_eshort+tcpu()-tt0
1703 #endif
1704 #endif
1705         call cartgrad
1706         call lagrangian
1707         if(.not.out1file .and. large) then
1708           write (iout,*) "energia_long",energia_long(0),
1709      &     " energia_short",energia_short(0),
1710      &     " total",energia_long(0)+energia_short(0)
1711           write (iout,*) "Initial fast-force accelerations"
1712           do i=0,nres
1713             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1714      &      (d_a(j,i+nres),j=1,3)
1715           enddo
1716         endif
1717 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1718         do i=0,2*nres
1719           do j=1,3
1720             d_a_short(j,i)=d_a(j,i)
1721           enddo
1722         enddo
1723 #ifdef MPI
1724         tt0=MPI_Wtime()
1725 #else
1726         tt0=tcpu()
1727 #endif
1728         call zerograd
1729         call etotal_long(energia_long)
1730         if (large) call enerprint(potEcomp)
1731 #ifdef TIMING_ENE
1732 #ifdef MPI
1733         t_elong=t_elong+MPI_Wtime()-tt0
1734 #else
1735         t_elong=t_elong+tcpu()-tt0
1736 #endif
1737 #endif
1738         call cartgrad
1739         call lagrangian
1740         if(.not.out1file .and. large) then
1741           write (iout,*) "energia_long",energia_long(0)
1742           write (iout,*) "Initial slow-force accelerations"
1743           do i=0,nres
1744             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1745      &      (d_a(j,i+nres),j=1,3)
1746           enddo
1747         endif
1748 #ifdef MPI
1749         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1750 #else
1751         t_enegrad=t_enegrad+tcpu()-tt0
1752 #endif
1753       endif
1754       return
1755       end
1756 c-----------------------------------------------------------
1757       subroutine random_vel
1758       implicit real*8 (a-h,o-z)
1759       include 'DIMENSIONS'
1760       include 'COMMON.CONTROL'
1761       include 'COMMON.VAR'
1762       include 'COMMON.MD'
1763 #ifndef LANG0
1764       include 'COMMON.LANGEVIN'
1765 #else
1766       include 'COMMON.LANGEVIN.lang0'
1767 #endif
1768       include 'COMMON.CHAIN'
1769       include 'COMMON.DERIV'
1770       include 'COMMON.GEO'
1771       include 'COMMON.LOCAL'
1772       include 'COMMON.INTERACT'
1773       include 'COMMON.IOUNITS'
1774       include 'COMMON.NAMES'
1775       include 'COMMON.TIME1'
1776       double precision xv,sigv,lowb,highb
1777 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1778 c First generate velocities in the eigenspace of the G matrix
1779 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1780 c      call flush(iout)
1781       xv=0.0d0
1782       ii=0
1783       do i=1,dimen
1784         do k=1,3
1785           ii=ii+1
1786           sigv=dsqrt((Rb*t_bath)/geigen(i))
1787           lowb=-5*sigv
1788           highb=5*sigv
1789           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1790 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1791 c     &      " d_t_work_new",d_t_work_new(ii)
1792         enddo
1793       enddo
1794 c diagnostics
1795 c      Ek1=0.0d0
1796 c      ii=0
1797 c      do i=1,dimen
1798 c        do k=1,3
1799 c          ii=ii+1
1800 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1801 c        enddo
1802 c      enddo
1803 c      write (iout,*) "Ek from eigenvectors",Ek1
1804 c end diagnostics
1805 c Transform velocities to UNRES coordinate space
1806       do k=0,2
1807         do i=1,dimen
1808           ind=(i-1)*3+k+1
1809           d_t_work(ind)=0.0d0
1810           do j=1,dimen
1811             d_t_work(ind)=d_t_work(ind)
1812      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1813           enddo
1814 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1815 c          call flush(iout)
1816         enddo
1817       enddo
1818 c Transfer to the d_t vector
1819       do j=1,3
1820         d_t(j,0)=d_t_work(j)
1821       enddo
1822       ind=3
1823       do i=nnt,nct-1
1824         do j=1,3
1825           ind=ind+1
1826           d_t(j,i)=d_t_work(ind)
1827         enddo
1828       enddo
1829       do i=nnt,nct
1830         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1831           do j=1,3
1832             ind=ind+1
1833             d_t(j,i+nres)=d_t_work(ind)
1834           enddo
1835         endif
1836       enddo
1837 c      call kinetic(EK)
1838 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1839 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1840 c      call flush(iout)
1841       return
1842       end
1843 #ifndef LANG0
1844 c-----------------------------------------------------------
1845       subroutine sd_verlet_p_setup
1846 c Sets up the parameters of stochastic Verlet algorithm
1847       implicit real*8 (a-h,o-z)
1848       include 'DIMENSIONS'
1849 #ifdef MPI
1850       include 'mpif.h'
1851 #endif
1852       include 'COMMON.CONTROL'
1853       include 'COMMON.VAR'
1854       include 'COMMON.MD'
1855 #ifndef LANG0
1856       include 'COMMON.LANGEVIN'
1857 #else
1858       include 'COMMON.LANGEVIN.lang0'
1859 #endif
1860       include 'COMMON.CHAIN'
1861       include 'COMMON.DERIV'
1862       include 'COMMON.GEO'
1863       include 'COMMON.LOCAL'
1864       include 'COMMON.INTERACT'
1865       include 'COMMON.IOUNITS'
1866       include 'COMMON.NAMES'
1867       include 'COMMON.TIME1'
1868       double precision emgdt(MAXRES6),
1869      & pterm,vterm,rho,rhoc,vsig,
1870      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1871      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1872      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1873       logical lprn /.false./
1874       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1875       double precision ktm
1876 #ifdef MPI
1877       tt0 = MPI_Wtime()
1878 #else
1879       tt0 = tcpu()
1880 #endif
1881 c
1882 c AL 8/17/04 Code adapted from tinker
1883 c
1884 c Get the frictional and random terms for stochastic dynamics in the
1885 c eigenspace of mass-scaled UNRES friction matrix
1886 c
1887       do i = 1, dimen
1888             gdt = fricgam(i) * d_time
1889 c
1890 c Stochastic dynamics reduces to simple MD for zero friction
1891 c
1892             if (gdt .le. zero) then
1893                pfric_vec(i) = 1.0d0
1894                vfric_vec(i) = d_time
1895                afric_vec(i) = 0.5d0 * d_time * d_time
1896                prand_vec(i) = 0.0d0
1897                vrand_vec1(i) = 0.0d0
1898                vrand_vec2(i) = 0.0d0
1899 c
1900 c Analytical expressions when friction coefficient is large
1901 c
1902             else
1903                if (gdt .ge. gdt_radius) then
1904                   egdt = dexp(-gdt)
1905                   pfric_vec(i) = egdt
1906                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1907                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1908                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1909                   vterm = 1.0d0 - egdt**2
1910                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1911 c
1912 c Use series expansions when friction coefficient is small
1913 c
1914                else
1915                   gdt2 = gdt * gdt
1916                   gdt3 = gdt * gdt2
1917                   gdt4 = gdt2 * gdt2
1918                   gdt5 = gdt2 * gdt3
1919                   gdt6 = gdt3 * gdt3
1920                   gdt7 = gdt3 * gdt4
1921                   gdt8 = gdt4 * gdt4
1922                   gdt9 = gdt4 * gdt5
1923                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
1924      &                          - gdt5/120.0d0 + gdt6/720.0d0
1925      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
1926      &                          - gdt9/362880.0d0) / fricgam(i)**2
1927                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
1928                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
1929                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
1930      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
1931      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
1932      &                       + 127.0d0*gdt9/90720.0d0
1933                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
1934      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
1935      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
1936      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
1937                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
1938      &                       - 17.0d0*gdt2/1280.0d0
1939      &                       + 17.0d0*gdt3/6144.0d0
1940      &                       + 40967.0d0*gdt4/34406400.0d0
1941      &                       - 57203.0d0*gdt5/275251200.0d0
1942      &                       - 1429487.0d0*gdt6/13212057600.0d0)
1943                end if
1944 c
1945 c Compute the scaling factors of random terms for the nonzero friction case
1946 c
1947                ktm = 0.5d0*d_time/fricgam(i)
1948                psig = dsqrt(ktm*pterm) / fricgam(i)
1949                vsig = dsqrt(ktm*vterm)
1950                rhoc = dsqrt(1.0d0 - rho*rho)
1951                prand_vec(i) = psig
1952                vrand_vec1(i) = vsig * rho
1953                vrand_vec2(i) = vsig * rhoc
1954             end if
1955       end do
1956       if (lprn) then
1957       write (iout,*)
1958      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
1959      &  " vrand_vec2"
1960       do i=1,dimen
1961         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
1962      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
1963       enddo
1964       endif
1965 c
1966 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
1967 c
1968 #ifndef   LANG0
1969       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
1970       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
1971       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
1972       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
1973       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
1974       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
1975 #endif
1976 #ifdef MPI
1977       t_sdsetup=t_sdsetup+MPI_Wtime()
1978 #else
1979       t_sdsetup=t_sdsetup+tcpu()-tt0
1980 #endif
1981       return
1982       end
1983 c-------------------------------------------------------------
1984       subroutine eigtransf1(n,ndim,ab,d,c)
1985       implicit none
1986       integer n,ndim
1987       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
1988       integer i,j,k
1989       do i=1,n
1990         do j=1,n
1991           c(i,j)=0.0d0
1992           do k=1,n
1993             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
1994           enddo
1995         enddo
1996       enddo
1997       return
1998       end
1999 c-------------------------------------------------------------
2000       subroutine eigtransf(n,ndim,a,b,d,c)
2001       implicit none
2002       integer n,ndim
2003       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2004       integer i,j,k
2005       do i=1,n
2006         do j=1,n
2007           c(i,j)=0.0d0
2008           do k=1,n
2009             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2010           enddo
2011         enddo
2012       enddo
2013       return
2014       end
2015 c-------------------------------------------------------------
2016       subroutine sd_verlet1
2017 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2018       implicit real*8 (a-h,o-z)
2019       include 'DIMENSIONS'
2020       include 'COMMON.CONTROL'
2021       include 'COMMON.VAR'
2022       include 'COMMON.MD'
2023 #ifndef LANG0
2024       include 'COMMON.LANGEVIN'
2025 #else
2026       include 'COMMON.LANGEVIN.lang0'
2027 #endif
2028       include 'COMMON.CHAIN'
2029       include 'COMMON.DERIV'
2030       include 'COMMON.GEO'
2031       include 'COMMON.LOCAL'
2032       include 'COMMON.INTERACT'
2033       include 'COMMON.IOUNITS'
2034       include 'COMMON.NAMES'
2035       double precision stochforcvec(MAXRES6)
2036       common /stochcalc/ stochforcvec
2037       logical lprn /.false./
2038
2039 c      write (iout,*) "dc_old"
2040 c      do i=0,nres
2041 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2042 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2043 c      enddo
2044       do j=1,3
2045         dc_work(j)=dc_old(j,0)
2046         d_t_work(j)=d_t_old(j,0)
2047         d_a_work(j)=d_a_old(j,0)
2048       enddo
2049       ind=3
2050       do i=nnt,nct-1
2051         do j=1,3
2052           dc_work(ind+j)=dc_old(j,i)
2053           d_t_work(ind+j)=d_t_old(j,i)
2054           d_a_work(ind+j)=d_a_old(j,i)
2055         enddo
2056         ind=ind+3
2057       enddo
2058       do i=nnt,nct
2059         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2060           do j=1,3
2061             dc_work(ind+j)=dc_old(j,i+nres)
2062             d_t_work(ind+j)=d_t_old(j,i+nres)
2063             d_a_work(ind+j)=d_a_old(j,i+nres)
2064           enddo
2065           ind=ind+3
2066         endif
2067       enddo
2068 #ifndef LANG0
2069       if (lprn) then
2070       write (iout,*)
2071      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2072      &  " vrand_mat2"
2073       do i=1,dimen
2074         do j=1,dimen
2075           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2076      &      vfric_mat(i,j),afric_mat(i,j),
2077      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2078         enddo
2079       enddo
2080       endif
2081       do i=1,dimen
2082         ddt1=0.0d0
2083         ddt2=0.0d0
2084         do j=1,dimen
2085           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2086      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2087           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2088           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2089         enddo
2090         d_t_work_new(i)=ddt1+0.5d0*ddt2
2091         d_t_work(i)=ddt1+ddt2
2092       enddo
2093 #endif
2094       do j=1,3
2095         dc(j,0)=dc_work(j)
2096         d_t(j,0)=d_t_work(j)
2097       enddo
2098       ind=3
2099       do i=nnt,nct-1
2100         do j=1,3
2101           dc(j,i)=dc_work(ind+j)
2102           d_t(j,i)=d_t_work(ind+j)
2103         enddo
2104         ind=ind+3
2105       enddo
2106       do i=nnt,nct
2107         if (itype(i).ne.10) then
2108           inres=i+nres
2109           do j=1,3
2110             dc(j,inres)=dc_work(ind+j)
2111             d_t(j,inres)=d_t_work(ind+j)
2112           enddo
2113           ind=ind+3
2114         endif
2115       enddo
2116       return
2117       end
2118 c--------------------------------------------------------------------------
2119       subroutine sd_verlet2
2120 c  Calculating the adjusted velocities for accelerations
2121       implicit real*8 (a-h,o-z)
2122       include 'DIMENSIONS'
2123       include 'COMMON.CONTROL'
2124       include 'COMMON.VAR'
2125       include 'COMMON.MD'
2126 #ifndef LANG0
2127       include 'COMMON.LANGEVIN'
2128 #else
2129       include 'COMMON.LANGEVIN.lang0'
2130 #endif
2131       include 'COMMON.CHAIN'
2132       include 'COMMON.DERIV'
2133       include 'COMMON.GEO'
2134       include 'COMMON.LOCAL'
2135       include 'COMMON.INTERACT'
2136       include 'COMMON.IOUNITS'
2137       include 'COMMON.NAMES'
2138       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2139       common /stochcalc/ stochforcvec
2140 c
2141 c Compute the stochastic forces which contribute to velocity change
2142 c
2143       call stochastic_force(stochforcvecV)
2144
2145 #ifndef LANG0
2146       do i=1,dimen
2147         ddt1=0.0d0
2148         ddt2=0.0d0
2149         do j=1,dimen
2150           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2151           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2152      &     vrand_mat2(i,j)*stochforcvecV(j)
2153         enddo
2154         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2155       enddo
2156 #endif
2157       do j=1,3
2158         d_t(j,0)=d_t_work(j)
2159       enddo
2160       ind=3
2161       do i=nnt,nct-1
2162         do j=1,3
2163           d_t(j,i)=d_t_work(ind+j)
2164         enddo
2165         ind=ind+3
2166       enddo
2167       do i=nnt,nct
2168         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2169           inres=i+nres
2170           do j=1,3
2171             d_t(j,inres)=d_t_work(ind+j)
2172           enddo
2173           ind=ind+3
2174         endif
2175       enddo
2176       return
2177       end
2178 c-----------------------------------------------------------
2179       subroutine sd_verlet_ciccotti_setup
2180 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2181 c version
2182       implicit real*8 (a-h,o-z)
2183       include 'DIMENSIONS'
2184 #ifdef MPI
2185       include 'mpif.h'
2186 #endif
2187       include 'COMMON.CONTROL'
2188       include 'COMMON.VAR'
2189       include 'COMMON.MD'
2190 #ifndef LANG0
2191       include 'COMMON.LANGEVIN'
2192 #else
2193       include 'COMMON.LANGEVIN.lang0'
2194 #endif
2195       include 'COMMON.CHAIN'
2196       include 'COMMON.DERIV'
2197       include 'COMMON.GEO'
2198       include 'COMMON.LOCAL'
2199       include 'COMMON.INTERACT'
2200       include 'COMMON.IOUNITS'
2201       include 'COMMON.NAMES'
2202       include 'COMMON.TIME1'
2203       double precision emgdt(MAXRES6),
2204      & pterm,vterm,rho,rhoc,vsig,
2205      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2206      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2207      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2208       logical lprn /.false./
2209       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2210       double precision ktm
2211 #ifdef MPI
2212       tt0 = MPI_Wtime()
2213 #else
2214       tt0 = tcpu()
2215 #endif
2216 c
2217 c AL 8/17/04 Code adapted from tinker
2218 c
2219 c Get the frictional and random terms for stochastic dynamics in the
2220 c eigenspace of mass-scaled UNRES friction matrix
2221 c
2222       do i = 1, dimen
2223             write (iout,*) "i",i," fricgam",fricgam(i)
2224             gdt = fricgam(i) * d_time
2225 c
2226 c Stochastic dynamics reduces to simple MD for zero friction
2227 c
2228             if (gdt .le. zero) then
2229                pfric_vec(i) = 1.0d0
2230                vfric_vec(i) = d_time
2231                afric_vec(i) = 0.5d0*d_time*d_time
2232                prand_vec(i) = afric_vec(i)
2233                vrand_vec2(i) = vfric_vec(i)
2234 c
2235 c Analytical expressions when friction coefficient is large
2236 c
2237             else
2238                egdt = dexp(-gdt)
2239                pfric_vec(i) = egdt
2240                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2241                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2242                prand_vec(i) = afric_vec(i)
2243                vrand_vec2(i) = vfric_vec(i)
2244 c
2245 c Compute the scaling factors of random terms for the nonzero friction case
2246 c
2247 c               ktm = 0.5d0*d_time/fricgam(i)
2248 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2249 c               vsig = dsqrt(ktm*vterm)
2250 c               prand_vec(i) = psig*afric_vec(i)
2251 c               vrand_vec2(i) = vsig*vfric_vec(i)
2252             end if
2253       end do
2254       if (lprn) then
2255       write (iout,*)
2256      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2257      &  " vrand_vec2"
2258       do i=1,dimen
2259         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2260      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2261       enddo
2262       endif
2263 c
2264 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2265 c
2266       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2267       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2268       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2269       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2270       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2271 #ifdef MPI
2272       t_sdsetup=t_sdsetup+MPI_Wtime()
2273 #else
2274       t_sdsetup=t_sdsetup+tcpu()-tt0
2275 #endif
2276       return
2277       end
2278 c-------------------------------------------------------------
2279       subroutine sd_verlet1_ciccotti
2280 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2281       implicit real*8 (a-h,o-z)
2282       include 'DIMENSIONS'
2283 #ifdef MPI
2284       include 'mpif.h'
2285 #endif
2286       include 'COMMON.CONTROL'
2287       include 'COMMON.VAR'
2288       include 'COMMON.MD'
2289 #ifndef LANG0
2290       include 'COMMON.LANGEVIN'
2291 #else
2292       include 'COMMON.LANGEVIN.lang0'
2293 #endif
2294       include 'COMMON.CHAIN'
2295       include 'COMMON.DERIV'
2296       include 'COMMON.GEO'
2297       include 'COMMON.LOCAL'
2298       include 'COMMON.INTERACT'
2299       include 'COMMON.IOUNITS'
2300       include 'COMMON.NAMES'
2301       double precision stochforcvec(MAXRES6)
2302       common /stochcalc/ stochforcvec
2303       logical lprn /.false./
2304
2305 c      write (iout,*) "dc_old"
2306 c      do i=0,nres
2307 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2308 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2309 c      enddo
2310       do j=1,3
2311         dc_work(j)=dc_old(j,0)
2312         d_t_work(j)=d_t_old(j,0)
2313         d_a_work(j)=d_a_old(j,0)
2314       enddo
2315       ind=3
2316       do i=nnt,nct-1
2317         do j=1,3
2318           dc_work(ind+j)=dc_old(j,i)
2319           d_t_work(ind+j)=d_t_old(j,i)
2320           d_a_work(ind+j)=d_a_old(j,i)
2321         enddo
2322         ind=ind+3
2323       enddo
2324       do i=nnt,nct
2325         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2326           do j=1,3
2327             dc_work(ind+j)=dc_old(j,i+nres)
2328             d_t_work(ind+j)=d_t_old(j,i+nres)
2329             d_a_work(ind+j)=d_a_old(j,i+nres)
2330           enddo
2331           ind=ind+3
2332         endif
2333       enddo
2334
2335 #ifndef LANG0
2336       if (lprn) then
2337       write (iout,*)
2338      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2339      &  " vrand_mat2"
2340       do i=1,dimen
2341         do j=1,dimen
2342           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2343      &      vfric_mat(i,j),afric_mat(i,j),
2344      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2345         enddo
2346       enddo
2347       endif
2348       do i=1,dimen
2349         ddt1=0.0d0
2350         ddt2=0.0d0
2351         do j=1,dimen
2352           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2353      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2354           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2355           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2356         enddo
2357         d_t_work_new(i)=ddt1+0.5d0*ddt2
2358         d_t_work(i)=ddt1+ddt2
2359       enddo
2360 #endif
2361       do j=1,3
2362         dc(j,0)=dc_work(j)
2363         d_t(j,0)=d_t_work(j)
2364       enddo
2365       ind=3
2366       do i=nnt,nct-1
2367         do j=1,3
2368           dc(j,i)=dc_work(ind+j)
2369           d_t(j,i)=d_t_work(ind+j)
2370         enddo
2371         ind=ind+3
2372       enddo
2373       do i=nnt,nct
2374         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2375           inres=i+nres
2376           do j=1,3
2377             dc(j,inres)=dc_work(ind+j)
2378             d_t(j,inres)=d_t_work(ind+j)
2379           enddo
2380           ind=ind+3
2381         endif
2382       enddo
2383       return
2384       end
2385 c--------------------------------------------------------------------------
2386       subroutine sd_verlet2_ciccotti
2387 c  Calculating the adjusted velocities for accelerations
2388       implicit real*8 (a-h,o-z)
2389       include 'DIMENSIONS'
2390       include 'COMMON.CONTROL'
2391       include 'COMMON.VAR'
2392       include 'COMMON.MD'
2393 #ifndef LANG0
2394       include 'COMMON.LANGEVIN'
2395 #else
2396       include 'COMMON.LANGEVIN.lang0'
2397 #endif
2398       include 'COMMON.CHAIN'
2399       include 'COMMON.DERIV'
2400       include 'COMMON.GEO'
2401       include 'COMMON.LOCAL'
2402       include 'COMMON.INTERACT'
2403       include 'COMMON.IOUNITS'
2404       include 'COMMON.NAMES'
2405       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2406       common /stochcalc/ stochforcvec
2407 c
2408 c Compute the stochastic forces which contribute to velocity change
2409 c
2410       call stochastic_force(stochforcvecV)
2411 #ifndef LANG0
2412       do i=1,dimen
2413         ddt1=0.0d0
2414         ddt2=0.0d0
2415         do j=1,dimen
2416
2417           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2418 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2419           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2420         enddo
2421         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2422       enddo
2423 #endif
2424       do j=1,3
2425         d_t(j,0)=d_t_work(j)
2426       enddo
2427       ind=3
2428       do i=nnt,nct-1
2429         do j=1,3
2430           d_t(j,i)=d_t_work(ind+j)
2431         enddo
2432         ind=ind+3
2433       enddo
2434       do i=nnt,nct
2435         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2436           inres=i+nres
2437           do j=1,3
2438             d_t(j,inres)=d_t_work(ind+j)
2439           enddo
2440           ind=ind+3
2441         endif
2442       enddo
2443       return
2444       end
2445 #endif