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