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