source/unres/src-HCD-5D/MD_A-MTS.F

   1       subroutine MD
   2 c------------------------------------------------
   3 c  The driver for molecular dynamics subroutines
   4 c------------------------------------------------
   5       implicit none
   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 #ifdef FIVEDIAG
  16        include 'COMMON.LAGRANGE.5diag'
  17 #else
  18        include 'COMMON.LAGRANGE'
  19 #endif
  20 #ifndef LANG0
  21       include 'COMMON.LANGEVIN'
  22 #else
  23 #ifdef FIVEDIAG
  24       include 'COMMON.LANGEVIN.lang0.5diag'
  25 #else
  26       include 'COMMON.LANGEVIN.lang0'
  27 #endif
  28 #endif
  29       include 'COMMON.CHAIN'
  30       include 'COMMON.DERIV'
  31       include 'COMMON.GEO'
  32       include 'COMMON.LOCAL'
  33       include 'COMMON.INTERACT'
  34       include 'COMMON.IOUNITS'
  35       include 'COMMON.NAMES'
  36       include 'COMMON.TIME1'
  37       include 'COMMON.HAIRPIN'
  38       double precision cm(3),L(3),vcm(3)
  39 #ifdef VOUT
  40       double precision v_work(maxres6),v_transf(maxres6)
  41 #endif
  42       integer ilen,rstcount
  43       external ilen
  44       character*50 tytul
  45       common /gucio/ cm
  46       integer itime
  47       logical ovrtim
  48       integer i,j,icount_scale,itime_scal
  49       integer nharp,iharp(4,maxres/3)
  50       double precision scalfac
  51       double precision tt0
  52 c
  53 #ifdef MPI
  54       if (ilen(tmpdir).gt.0)
  55      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
  56      &        //liczba(:ilen(liczba))//'.rst')
  57 #else
  58       if (ilen(tmpdir).gt.0)
  59      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
  60 #endif
  61       write (iout,*) "MD lang",lang
  62       t_MDsetup=0.0d0
  63       t_langsetup=0.0d0
  64       t_MD=0.0d0
  65       t_enegrad=0.0d0
  66       t_sdsetup=0.0d0
  67       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  68 #ifdef MPI
  69       tt0=MPI_Wtime()
  70 #else
  71       tt0 = tcpu()
  72 #endif
  73 c Determine the inverse of the inertia matrix.
  74       call setup_MD_matrices
  75 c Initialize MD
  76       call init_MD
  77 #ifdef MPI
  78       t_MDsetup = MPI_Wtime()-tt0
  79 #else
  80       t_MDsetup = tcpu()-tt0
  81 #endif
  82       rstcount=0
  83 c   Entering the MD loop
  84 #ifdef MPI
  85       tt0 = MPI_Wtime()
  86 #else
  87       tt0 = tcpu()
  88 #endif
  89       if (lang.eq.2 .or. lang.eq.3) then
  90 #ifndef   LANG0
  91         call setup_fricmat
  92         if (lang.eq.2) then
  93           call sd_verlet_p_setup
  94         else
  95           call sd_verlet_ciccotti_setup
  96         endif
  97         do i=1,dimen
  98           do j=1,dimen
  99             pfric0_mat(i,j,0)=pfric_mat(i,j)
 100             afric0_mat(i,j,0)=afric_mat(i,j)
 101             vfric0_mat(i,j,0)=vfric_mat(i,j)
 102             prand0_mat(i,j,0)=prand_mat(i,j)
 103             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 104             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 105           enddo
 106         enddo
 107         flag_stoch(0)=.true.
 108         do i=1,maxflag_stoch
 109           flag_stoch(i)=.false.
 110         enddo
 111 #else
 112         write (iout,*)
 113      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 114 #ifdef MPI
 115         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 116 #endif
 117         stop
 118 #endif
 119       else if (lang.eq.1 .or. lang.eq.4) then
 120         call setup_fricmat
 121       endif
 122 #ifdef MPI
 123       t_langsetup=MPI_Wtime()-tt0
 124       tt0=MPI_Wtime()
 125 #else
 126       t_langsetup=tcpu()-tt0
 127       tt0=tcpu()
 128 #endif
 129       do itime=1,n_timestep
 130         if (ovrtim()) exit
 131         if (large.and. mod(itime,ntwe).eq.0)
 132      &    write (iout,*) "itime",itime
 133         rstcount=rstcount+1
 134         if (lang.gt.0 .and. surfarea .and.
 135      &      mod(itime,reset_fricmat).eq.0) then
 136           if (lang.eq.2 .or. lang.eq.3) then
 137 #ifndef LANG0
 138             call setup_fricmat
 139             if (lang.eq.2) then
 140               call sd_verlet_p_setup
 141             else
 142               call sd_verlet_ciccotti_setup
 143             endif
 144             do i=1,dimen
 145               do j=1,dimen
 146                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 147                 afric0_mat(i,j,0)=afric_mat(i,j)
 148                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 149                 prand0_mat(i,j,0)=prand_mat(i,j)
 150                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 151                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 152               enddo
 153             enddo
 154             flag_stoch(0)=.true.
 155             do i=1,maxflag_stoch
 156               flag_stoch(i)=.false.
 157             enddo
 158 #endif
 159           else if (lang.eq.1 .or. lang.eq.4) then
 160             call setup_fricmat
 161           endif
 162           write (iout,'(a,i10)')
 163      &      "Friction matrix reset based on surface area, itime",itime
 164         endif
 165         if (reset_vel .and. tbf .and. lang.eq.0
 166      &      .and. mod(itime,count_reset_vel).eq.0) then
 167           call random_vel
 168           write(iout,'(a,f20.2)')
 169      &     "Velocities reset to random values, time",totT
 170           do i=0,2*nres
 171             do j=1,3
 172               d_t_old(j,i)=d_t(j,i)
 173             enddo
 174           enddo
 175         endif
 176         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 177           call inertia_tensor
 178           call vcm_vel(vcm)
 179           do j=1,3
 180              d_t(j,0)=d_t(j,0)-vcm(j)
 181           enddo
 182           call kinetic(EK)
 183           kinetic_T=2.0d0/(dimen3*Rb)*EK
 184           scalfac=dsqrt(T_bath/kinetic_T)
 185           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 186           do i=0,2*nres
 187             do j=1,3
 188               d_t_old(j,i)=scalfac*d_t(j,i)
 189             enddo
 190           enddo
 191         endif
 192         if (lang.ne.4) then
 193           if (RESPA) then
 194 c Time-reversible RESPA algorithm
 195 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 196             call RESPA_step(itime)
 197           else
 198 c Variable time step algorithm.
 199             call velverlet_step(itime)
 200           endif
 201         else
 202 #ifdef BROWN
 203           call brown_step(itime)
 204 #else
 205           print *,"Brown dynamics not here!"
 206 #ifdef MPI
 207           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 208 #endif
 209           stop
 210 #endif
 211         endif
 212         if (ntwe.ne.0) then
 213          if (mod(itime,ntwe).eq.0) then
 214            call statout(itime)
 215 C           call enerprint(potEcomp)
 216 C           print *,itime,'AFM',Eafmforc,etot
 217          endif
 218 #ifdef VOUT
 219         do j=1,3
 220           v_work(j)=d_t(j,0)
 221         enddo
 222         ind=3
 223         do i=nnt,nct-1
 224           do j=1,3
 225             ind=ind+1
 226             v_work(ind)=d_t(j,i)
 227           enddo
 228         enddo
 229         do i=nnt,nct
 230           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 231             do j=1,3
 232               ind=ind+1
 233               v_work(ind)=d_t(j,i+nres)
 234             enddo
 235           endif
 236         enddo
 237
 238         write (66,'(80f10.5)')
 239      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 240         do i=1,ind
 241           v_transf(i)=0.0d0
 242           do j=1,ind
 243             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 244           enddo
 245            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 246         enddo
 247         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 248 #endif
 249         endif
 250         if (mod(itime,ntwx).eq.0) then
 251 c          write(iout,*) 'time=',itime
 252 C          call check_ecartint
 253           call returnbox
 254           write (tytul,'("time",f8.2)') totT
 255           if(mdpdb) then
 256              call hairpin(.true.,nharp,iharp)
 257              call secondary2(.true.)
 258              call pdbout(potE,tytul,ipdb)
 259           else
 260              call cartout(totT)
 261           endif
 262         endif
 263         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
 264            open(irest2,file=rest2name,status='unknown')
 265            write(irest2,*) totT,EK,potE,totE,t_bath
 266            do i=1,2*nres
 267             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 268            enddo
 269            do i=1,2*nres
 270             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 271            enddo
 272           close(irest2)
 273           rstcount=0
 274         endif
 275       enddo
 276 #ifdef MPI
 277       t_MD=MPI_Wtime()-tt0
 278 #else
 279       t_MD=tcpu()-tt0
 280 #endif
 281       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 282      &  '  Timing  ',
 283      & 'MD calculations setup:',t_MDsetup,
 284      & 'Energy & gradient evaluation:',t_enegrad,
 285      & 'Stochastic MD setup:',t_langsetup,
 286      & 'Stochastic MD step setup:',t_sdsetup,
 287      & 'MD steps:',t_MD
 288       write (iout,'(/28(1h=),a25,27(1h=))')
 289      & '  End of MD calculation  '
 290 #ifdef TIMING_ENE
 291       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
 292      &  " eshort",t_eshort
 293       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
 294      & " time_fricmatmult",time_fricmatmult," time_fsample ",
 295      & time_fsample
 296 #endif
 297       return
 298       end
 299 c-------------------------------------------------------------------------------
 300       subroutine velverlet_step(itime)
 301 c-------------------------------------------------------------------------------
 302 c  Perform a single velocity Verlet step; the time step can be rescaled if
 303 c  increments in accelerations exceed the threshold
 304 c-------------------------------------------------------------------------------
 305       implicit none
 306       include 'DIMENSIONS'
 307 #ifdef MPI
 308       include 'mpif.h'
 309       integer ierror,ierrcode,errcode
 310 #endif
 311       include 'COMMON.SETUP'
 312       include 'COMMON.CONTROL'
 313       include 'COMMON.VAR'
 314       include 'COMMON.MD'
 315 #ifdef FIVEDIAG
 316        include 'COMMON.LAGRANGE.5diag'
 317 #else
 318        include 'COMMON.LAGRANGE'
 319 #endif
 320 #ifndef LANG0
 321       include 'COMMON.LANGEVIN'
 322 #else
 323 #ifdef FIVEDIAG
 324       include 'COMMON.LANGEVIN.lang0.5diag'
 325 #else
 326       include 'COMMON.LANGEVIN.lang0'
 327 #endif
 328 #endif
 329       include 'COMMON.CHAIN'
 330       include 'COMMON.DERIV'
 331       include 'COMMON.GEO'
 332       include 'COMMON.LOCAL'
 333       include 'COMMON.INTERACT'
 334       include 'COMMON.IOUNITS'
 335       include 'COMMON.NAMES'
 336       include 'COMMON.TIME1'
 337       include 'COMMON.MUCA'
 338       double precision vcm(3),incr(3)
 339       double precision cm(3),L(3)
 340       integer ilen,count,rstcount
 341       external ilen
 342       character*50 tytul
 343       integer maxcount_scale /20/
 344       common /gucio/ cm
 345       double precision stochforcvec(MAXRES6)
 346       common /stochcalc/ stochforcvec
 347       integer itime,icount_scale,itime_scal,ifac_time,i,j,itt
 348       logical scale
 349       double precision epdrift,fac_time
 350       double precision tt0
 351 c
 352       scale=.true.
 353       icount_scale=0
 354       if (lang.eq.1) then
 355         call sddir_precalc
 356       else if (lang.eq.2 .or. lang.eq.3) then
 357 #ifndef LANG0
 358         call stochastic_force(stochforcvec)
 359 #else
 360         write (iout,*)
 361      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 362 #ifdef MPI
 363         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 364 #endif
 365         stop
 366 #endif
 367       endif
 368       itime_scal=0
 369       do while (scale)
 370         icount_scale=icount_scale+1
 371         if (icount_scale.gt.maxcount_scale) then
 372           write (iout,*)
 373      &      "ERROR: too many attempts at scaling down the time step. ",
 374      &      "amax=",amax,"epdrift=",epdrift,
 375      &      "damax=",damax,"edriftmax=",edriftmax,
 376      &      "d_time=",d_time
 377           call flush(iout)
 378 #ifdef MPI
 379           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
 380 #endif
 381           stop
 382         endif
 383 c First step of the velocity Verlet algorithm
 384         if (lang.eq.2) then
 385 #ifndef LANG0
 386           call sd_verlet1
 387 #endif
 388         else if (lang.eq.3) then
 389 #ifndef LANG0
 390           call sd_verlet1_ciccotti
 391 #endif
 392         else if (lang.eq.1) then
 393           call sddir_verlet1
 394         else
 395           call verlet1
 396         endif
 397 c Build the chain from the newly calculated coordinates
 398         call chainbuild_cart
 399         if (rattle) call rattle1
 400         if (ntwe.ne.0) then
 401         if (large.and. mod(itime,ntwe).eq.0) then
 402           write (iout,*) "Cartesian and internal coordinates: step 1"
 403           call cartprint
 404           call intout
 405           write (iout,*) "dC"
 406           do i=0,nres
 407             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
 408      &      (dc(j,i+nres),j=1,3)
 409           enddo
 410           write (iout,*) "Accelerations"
 411           do i=0,nres
 412             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 413      &      (d_a(j,i+nres),j=1,3)
 414           enddo
 415           write (iout,*) "Velocities, step 1"
 416           do i=0,nres
 417             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 418      &      (d_t(j,i+nres),j=1,3)
 419           enddo
 420         endif
 421         endif
 422 #ifdef MPI
 423         tt0 = MPI_Wtime()
 424 #else
 425         tt0 = tcpu()
 426 #endif
 427 c Calculate energy and forces
 428         call zerograd
 429         call etotal(potEcomp)
 430 ! AL 4/17/17: Reduce the steps if NaNs occurred.
 431         if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
 432           d_time=d_time/2
 433           cycle
 434         endif
 435 ! end change
 436         if (large.and. mod(itime,ntwe).eq.0)
 437      &    call enerprint(potEcomp)
 438 #ifdef TIMING_ENE
 439 #ifdef MPI
 440         t_etotal=t_etotal+MPI_Wtime()-tt0
 441 #else
 442         t_etotal=t_etotal+tcpu()-tt0
 443 #endif
 444 #endif
 445         potE=potEcomp(0)-potEcomp(27)
 446         call cartgrad
 447 c Get the new accelerations
 448         call lagrangian
 449 #ifdef MPI
 450         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 451 #else
 452         t_enegrad=t_enegrad+tcpu()-tt0
 453 #endif
 454 c Determine maximum acceleration and scale down the timestep if needed
 455         call max_accel
 456         amax=amax/(itime_scal+1)**2
 457         call predict_edrift(epdrift)
 458         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
 459 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 460           scale=.true.
 461           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 462      &      /dlog(2.0d0)+1
 463           itime_scal=itime_scal+ifac_time
 464 c          fac_time=dmin1(damax/amax,0.5d0)
 465           fac_time=0.5d0**ifac_time
 466           d_time=d_time*fac_time
 467           if (lang.eq.2 .or. lang.eq.3) then
 468 #ifndef LANG0
 469 c            write (iout,*) "Calling sd_verlet_setup: 1"
 470 c Rescale the stochastic forces and recalculate or restore
 471 c the matrices of tinker integrator
 472             if (itime_scal.gt.maxflag_stoch) then
 473               if (large) write (iout,'(a,i5,a)')
 474      &         "Calculate matrices for stochastic step;",
 475      &         " itime_scal ",itime_scal
 476               if (lang.eq.2) then
 477                 call sd_verlet_p_setup
 478               else
 479                 call sd_verlet_ciccotti_setup
 480               endif
 481               write (iout,'(2a,i3,a,i3,1h.)')
 482      &         "Warning: cannot store matrices for stochastic",
 483      &         " integration because the index",itime_scal,
 484      &         " is greater than",maxflag_stoch
 485               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 486      &         " integration Langevin algorithm for better efficiency."
 487             else if (flag_stoch(itime_scal)) then
 488               if (large) write (iout,'(a,i5,a,l1)')
 489      &         "Restore matrices for stochastic step; itime_scal ",
 490      &         itime_scal," flag ",flag_stoch(itime_scal)
 491               do i=1,dimen
 492                 do j=1,dimen
 493                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 494                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 495                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 496                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 497                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 498                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 499                 enddo
 500               enddo
 501             else
 502               if (large) write (iout,'(2a,i5,a,l1)')
 503      &         "Calculate & store matrices for stochastic step;",
 504      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 505               if (lang.eq.2) then
 506                 call sd_verlet_p_setup
 507               else
 508                 call sd_verlet_ciccotti_setup
 509               endif
 510               flag_stoch(ifac_time)=.true.
 511               do i=1,dimen
 512                 do j=1,dimen
 513                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 514                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 515                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 516                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 517                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 518                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 519                 enddo
 520               enddo
 521             endif
 522             fac_time=1.0d0/dsqrt(fac_time)
 523             do i=1,dimen
 524               stochforcvec(i)=fac_time*stochforcvec(i)
 525             enddo
 526 #endif
 527           else if (lang.eq.1) then
 528 c Rescale the accelerations due to stochastic forces
 529             fac_time=1.0d0/dsqrt(fac_time)
 530             do i=1,dimen
 531               d_as_work(i)=d_as_work(i)*fac_time
 532             enddo
 533           endif
 534           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 535      &      "itime",itime," Timestep scaled down to ",
 536      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 537         else
 538 c Second step of the velocity Verlet algorithm
 539           if (lang.eq.2) then
 540 #ifndef LANG0
 541             call sd_verlet2
 542 #endif
 543           else if (lang.eq.3) then
 544 #ifndef LANG0
 545             call sd_verlet2_ciccotti
 546 #endif
 547           else if (lang.eq.1) then
 548             call sddir_verlet2
 549           else
 550             call verlet2
 551           endif
 552           if (rattle) call rattle2
 553           totT=totT+d_time
 554           totTafm=totT
 555 C          print *,totTafm,"TU?"
 556           if (d_time.ne.d_time0) then
 557             d_time=d_time0
 558 #ifndef   LANG0
 559             if (lang.eq.2 .or. lang.eq.3) then
 560               if (large) write (iout,'(a)')
 561      &         "Restore original matrices for stochastic step"
 562 c              write (iout,*) "Calling sd_verlet_setup: 2"
 563 c Restore the matrices of tinker integrator if the time step has been restored
 564               do i=1,dimen
 565                 do j=1,dimen
 566                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 567                   afric_mat(i,j)=afric0_mat(i,j,0)
 568                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 569                   prand_mat(i,j)=prand0_mat(i,j,0)
 570                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 571                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 572                 enddo
 573               enddo
 574             endif
 575 #endif
 576           endif
 577           scale=.false.
 578         endif
 579       enddo
 580 c Calculate the kinetic and the total energy and the kinetic temperature
 581       call kinetic(EK)
 582       totE=EK+potE
 583 c diagnostics
 584 c      call kinetic1(EK1)
 585 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 586 c end diagnostics
 587 c Couple the system to Berendsen bath if needed
 588       if (tbf .and. lang.eq.0) then
 589         call verlet_bath
 590       endif
 591       kinetic_T=2.0d0/(dimen3*Rb)*EK
 592 c Backup the coordinates, velocities, and accelerations
 593       do i=0,2*nres
 594         do j=1,3
 595           dc_old(j,i)=dc(j,i)
 596           d_t_old(j,i)=d_t(j,i)
 597           d_a_old(j,i)=d_a(j,i)
 598         enddo
 599       enddo
 600       if (ntwe.ne.0) then
 601       if (mod(itime,ntwe).eq.0 .and. large) then
 602         write (iout,*) "Velocities, step 2"
 603         do i=0,nres
 604           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 605      &      (d_t(j,i+nres),j=1,3)
 606         enddo
 607       endif
 608       endif
 609       return
 610       end
 611 c-------------------------------------------------------------------------------
 612       subroutine RESPA_step(itime)
 613 c-------------------------------------------------------------------------------
 614 c  Perform a single RESPA step.
 615 c-------------------------------------------------------------------------------
 616       implicit none
 617       include 'DIMENSIONS'
 618 #ifdef MPI
 619       include 'mpif.h'
 620       integer IERROR,ERRCODE
 621 #endif
 622       include 'COMMON.SETUP'
 623       include 'COMMON.CONTROL'
 624       include 'COMMON.VAR'
 625       include 'COMMON.MD'
 626 #ifdef FIVEDIAG
 627        include 'COMMON.LAGRANGE.5diag'
 628 #else
 629        include 'COMMON.LAGRANGE'
 630 #endif
 631 #ifndef LANG0
 632       include 'COMMON.LANGEVIN'
 633 #else
 634 #ifdef FIVEDIAG
 635       include 'COMMON.LANGEVIN.lang0.5diag'
 636 #else
 637       include 'COMMON.LANGEVIN.lang0'
 638 #endif
 639 #endif
 640       include 'COMMON.CHAIN'
 641       include 'COMMON.DERIV'
 642       include 'COMMON.GEO'
 643       include 'COMMON.LOCAL'
 644       include 'COMMON.INTERACT'
 645       include 'COMMON.IOUNITS'
 646       include 'COMMON.NAMES'
 647       include 'COMMON.TIME1'
 648       logical lprint_short
 649       common /shortcheck/ lprint_short
 650       double precision energia_short(0:n_ene),
 651      & energia_long(0:n_ene)
 652       double precision cm(3),L(3),vcm(3),incr(3)
 653       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
 654      & d_a_old0(3,0:maxres2)
 655       integer i,j
 656       double precision fac_time
 657       logical PRINT_AMTS_MSG /.false./
 658       integer ilen,count,rstcount
 659       external ilen
 660       character*50 tytul
 661       integer maxcount_scale /10/
 662       common /gucio/ cm
 663       double precision stochforcvec(MAXRES6)
 664       common /stochcalc/ stochforcvec
 665       integer itime
 666       logical scale
 667       integer itt
 668       common /cipiszcze/ itt
 669       integer itsplit
 670       double precision epdrift,epdriftmax
 671       double precision tt0
 672       itt=itime
 673       if (ntwe.ne.0) then
 674       if (large.and. mod(itime,ntwe).eq.0) then
 675         write (iout,*) "***************** RESPA itime",itime
 676         write (iout,*) "Cartesian and internal coordinates: step 0"
 677 c        call cartprint
 678         call pdbout(0.0d0,
 679      &  "cipiszcze                                         ",iout)
 680         call intout
 681         write (iout,*) "Accelerations from long-range forces"
 682         do i=0,nres
 683           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 684      &      (d_a(j,i+nres),j=1,3)
 685         enddo
 686         write (iout,*) "Velocities, step 0"
 687         do i=0,nres
 688           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 689      &      (d_t(j,i+nres),j=1,3)
 690         enddo
 691       endif
 692       endif
 693 c
 694 c Perform the initial RESPA step (increment velocities)
 695 c      write (iout,*) "*********************** RESPA ini"
 696       call RESPA_vel
 697       if (ntwe.ne.0) then
 698       if (mod(itime,ntwe).eq.0 .and. large) then
 699         write (iout,*) "Velocities, end"
 700         do i=0,nres
 701           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 702      &      (d_t(j,i+nres),j=1,3)
 703         enddo
 704       endif
 705       endif
 706 c Compute the short-range forces
 707 #ifdef MPI
 708       tt0 =MPI_Wtime()
 709 #else
 710       tt0 = tcpu()
 711 #endif
 712 C 7/2/2009 commented out
 713 c      call zerograd
 714 c      call etotal_short(energia_short)
 715 c      call cartgrad
 716 c      call lagrangian
 717 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
 718         do i=0,2*nres
 719           do j=1,3
 720             d_a(j,i)=d_a_short(j,i)
 721           enddo
 722         enddo
 723       if (ntwe.ne.0) then
 724       if (large.and. mod(itime,ntwe).eq.0) then
 725         write (iout,*) "energia_short",energia_short(0)
 726         write (iout,*) "Accelerations from short-range forces"
 727         do i=0,nres
 728           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 729      &      (d_a(j,i+nres),j=1,3)
 730         enddo
 731       endif
 732       endif
 733 #ifdef MPI
 734         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 735 #else
 736         t_enegrad=t_enegrad+tcpu()-tt0
 737 #endif
 738       do i=0,2*nres
 739         do j=1,3
 740           dc_old(j,i)=dc(j,i)
 741           d_t_old(j,i)=d_t(j,i)
 742           d_a_old(j,i)=d_a(j,i)
 743         enddo
 744       enddo
 745 c 6/30/08 A-MTS: attempt at increasing the split number
 746       do i=0,2*nres
 747         do j=1,3
 748           dc_old0(j,i)=dc_old(j,i)
 749           d_t_old0(j,i)=d_t_old(j,i)
 750           d_a_old0(j,i)=d_a_old(j,i)
 751         enddo
 752       enddo
 753       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 754       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 755 c
 756       scale=.true.
 757       d_time0=d_time
 758       do while (scale)
 759
 760       scale=.false.
 761 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 762 c Split the time step
 763       d_time=d_time0/ntime_split
 764 c Perform the short-range RESPA steps (velocity Verlet increments of
 765 c positions and velocities using short-range forces)
 766 c      write (iout,*) "*********************** RESPA split"
 767       do itsplit=1,ntime_split
 768         if (lang.eq.1) then
 769           call sddir_precalc
 770         else if (lang.eq.2 .or. lang.eq.3) then
 771 #ifndef LANG0
 772           call stochastic_force(stochforcvec)
 773 #else
 774           write (iout,*)
 775      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 776 #ifdef MPI
 777           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 778 #endif
 779           stop
 780 #endif
 781         endif
 782 c First step of the velocity Verlet algorithm
 783         if (lang.eq.2) then
 784 #ifndef LANG0
 785           call sd_verlet1
 786 #endif
 787         else if (lang.eq.3) then
 788 #ifndef LANG0
 789           call sd_verlet1_ciccotti
 790 #endif
 791         else if (lang.eq.1) then
 792           call sddir_verlet1
 793         else
 794           call verlet1
 795         endif
 796 c Build the chain from the newly calculated coordinates
 797         call chainbuild_cart
 798         if (rattle) call rattle1
 799         if (ntwe.ne.0) then
 800         if (large.and. mod(itime,ntwe).eq.0) then
 801           write (iout,*) "***** ITSPLIT",itsplit
 802           write (iout,*) "Cartesian and internal coordinates: step 1"
 803           call pdbout(0.0d0,
 804      &     "cipiszcze                                         ",iout)
 805 c          call cartprint
 806           call intout
 807           write (iout,*) "Velocities, step 1"
 808           do i=0,nres
 809             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 810      &        (d_t(j,i+nres),j=1,3)
 811           enddo
 812         endif
 813         endif
 814 #ifdef MPI
 815         tt0 = MPI_Wtime()
 816 #else
 817         tt0 = tcpu()
 818 #endif
 819 c Calculate energy and forces
 820 c        if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
 821         call zerograd
 822         call etotal_short(energia_short)
 823 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
 824         if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) )
 825      &  then
 826           if (PRINT_AMTS_MSG) write (iout,*)
 827      &     "Infinities/NaNs in energia_short",
 828      &      energia_short(0),"; increasing ntime_split to",ntime_split
 829           ntime_split=ntime_split*2
 830           if (ntime_split.gt.maxtime_split) then
 831 #ifdef MPI
 832           write (iout,*) "Cannot rescue the run; aborting job.",
 833      &      " Retry with a smaller time step"
 834           call flush(iout)
 835           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 836 #else
 837           write (iout,*) "Cannot rescue the run; terminating.",
 838      &      " Retry with a smaller time step"
 839 #endif
 840           endif
 841           exit
 842         endif
 843 ! End change
 844         if (large.and. mod(itime,ntwe).eq.0)
 845      &    call enerprint(energia_short)
 846 #ifdef TIMING_ENE
 847 #ifdef MPI
 848         t_eshort=t_eshort+MPI_Wtime()-tt0
 849 #else
 850         t_eshort=t_eshort+tcpu()-tt0
 851 #endif
 852 #endif
 853         call cartgrad
 854 c Get the new accelerations
 855         call lagrangian
 856 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 857         do i=0,2*nres
 858           do j=1,3
 859             d_a_short(j,i)=d_a(j,i)
 860           enddo
 861         enddo
 862         if (ntwe.ne.0) then
 863         if (large.and. mod(itime,ntwe).eq.0) then
 864           write (iout,*)"energia_short",energia_short(0)
 865           write (iout,*) "Accelerations from short-range forces"
 866           do i=0,nres
 867             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 868      &        (d_a(j,i+nres),j=1,3)
 869           enddo
 870         endif
 871         endif
 872 c 6/30/08 A-MTS
 873 c Determine maximum acceleration and scale down the timestep if needed
 874         call max_accel
 875         amax=amax/ntime_split**2
 876         call predict_edrift(epdrift)
 877         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 878      &   write (iout,*) "amax",amax," damax",damax,
 879      &   " epdrift",epdrift," epdriftmax",epdriftmax
 880 c Exit loop and try with increased split number if the change of
 881 c acceleration is too big
 882         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 883           if (ntime_split.lt.maxtime_split) then
 884             scale=.true.
 885             ntime_split=ntime_split*2
 886             do i=0,2*nres
 887               do j=1,3
 888                 dc_old(j,i)=dc_old0(j,i)
 889                 d_t_old(j,i)=d_t_old0(j,i)
 890                 d_a_old(j,i)=d_a_old0(j,i)
 891               enddo
 892             enddo
 893             if (PRINT_AMTS_MSG) then
 894             write (iout,*) "acceleration/energy drift too large",amax,
 895      &      epdrift," split increased to ",ntime_split," itime",itime,
 896      &       " itsplit",itsplit
 897             endif
 898             exit
 899           else
 900             write (iout,*)
 901      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 902      &       maxtime_split,
 903      &      " already reached!!! Trying to carry on!"
 904           endif
 905         endif
 906 #ifdef MPI
 907         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 908 #else
 909         t_enegrad=t_enegrad+tcpu()-tt0
 910 #endif
 911 c Second step of the velocity Verlet algorithm
 912         if (lang.eq.2) then
 913 #ifndef LANG0
 914           call sd_verlet2
 915 #endif
 916         else if (lang.eq.3) then
 917 #ifndef LANG0
 918           call sd_verlet2_ciccotti
 919 #endif
 920         else if (lang.eq.1) then
 921           call sddir_verlet2
 922         else
 923           call verlet2
 924         endif
 925         if (rattle) call rattle2
 926 c Backup the coordinates, velocities, and accelerations
 927         do i=0,2*nres
 928           do j=1,3
 929             dc_old(j,i)=dc(j,i)
 930             d_t_old(j,i)=d_t(j,i)
 931             d_a_old(j,i)=d_a(j,i)
 932           enddo
 933         enddo
 934       enddo
 935
 936       enddo ! while scale
 937
 938 c Restore the time step
 939       d_time=d_time0
 940 c Compute long-range forces
 941 #ifdef MPI
 942       tt0 =MPI_Wtime()
 943 #else
 944       tt0 = tcpu()
 945 #endif
 946       call zerograd
 947       call etotal_long(energia_long)
 948 ! AL 4/17/2017 Handling NaNs
 949       if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
 950 #ifdef MPI
 951         write (iout,*)
 952      &   "Infinitied/NaNs in energia_long, Aborting MPI job."
 953         call enerprint(energia_long)
 954         call flush(iout)
 955         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 956 #else
 957         write (iout,*) "Infinitied/NaNs in energia_long, terminating."
 958         call enerprint(energia_long)
 959         stop
 960 #endif
 961       endif
 962 ! end change
 963 c        lprint_short=.false.
 964       if (large.and. mod(itime,ntwe).eq.0)
 965      &    call enerprint(energia_long)
 966 #ifdef TIMING_ENE
 967 #ifdef MPI
 968         t_elong=t_elong+MPI_Wtime()-tt0
 969 #else
 970         t_elong=t_elong+tcpu()-tt0
 971 #endif
 972 #endif
 973       call cartgrad
 974       call lagrangian
 975 #ifdef MPI
 976         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 977 #else
 978         t_enegrad=t_enegrad+tcpu()-tt0
 979 #endif
 980 c Compute accelerations from long-range forces
 981       if (ntwe.ne.0) then
 982       if (large.and. mod(itime,ntwe).eq.0) then
 983         write (iout,*) "energia_long",energia_long(0)
 984         write (iout,*) "Cartesian and internal coordinates: step 2"
 985 c        call cartprint
 986         call pdbout(0.0d0,
 987      &  'cipiszcze                                         ',iout)
 988         call intout
 989         write (iout,*) "Accelerations from long-range forces"
 990         do i=0,nres
 991           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 992      &      (d_a(j,i+nres),j=1,3)
 993         enddo
 994         write (iout,*) "Velocities, step 2"
 995         do i=0,nres
 996           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 997      &      (d_t(j,i+nres),j=1,3)
 998         enddo
 999       endif
1000       endif
1001 c Compute the final RESPA step (increment velocities)
1002 c      write (iout,*) "*********************** RESPA fin"
1003       call RESPA_vel
1004 c Compute the complete potential energy
1005       do i=0,n_ene
1006         potEcomp(i)=energia_short(i)+energia_long(i)
1007       enddo
1008       potE=potEcomp(0)-potEcomp(27)
1009       if (ntwe.ne.0) then
1010       if (large.and. mod(itime,ntwe).eq.0) then
1011         call enerprint(potEcomp)
1012         write (iout,*) "potE",potE
1013       endif
1014       endif
1015 c      potE=energia_short(0)+energia_long(0)
1016       totT=totT+d_time
1017       totTafm=totT
1018 c Calculate the kinetic and the total energy and the kinetic temperature
1019       call kinetic(EK)
1020       totE=EK+potE
1021 c Couple the system to Berendsen bath if needed
1022       if (tbf .and. lang.eq.0) then
1023         call verlet_bath
1024       endif
1025       kinetic_T=2.0d0/(dimen3*Rb)*EK
1026 c Backup the coordinates, velocities, and accelerations
1027       if (ntwe.ne.0) then
1028       if (mod(itime,ntwe).eq.0 .and. large) then
1029         write (iout,*) "Velocities, end"
1030         do i=0,nres
1031           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1032      &      (d_t(j,i+nres),j=1,3)
1033         enddo
1034       endif
1035       endif
1036       return
1037       end
1038 c---------------------------------------------------------------------
1039       subroutine RESPA_vel
1040 c  First and last RESPA step (incrementing velocities using long-range
1041 c  forces).
1042       implicit none
1043       include 'DIMENSIONS'
1044       include 'COMMON.CONTROL'
1045       include 'COMMON.VAR'
1046       include 'COMMON.MD'
1047 #ifdef FIVEDIAG
1048        include 'COMMON.LAGRANGE.5diag'
1049 #else
1050        include 'COMMON.LAGRANGE'
1051 #endif
1052       include 'COMMON.CHAIN'
1053       include 'COMMON.DERIV'
1054       include 'COMMON.GEO'
1055       include 'COMMON.LOCAL'
1056       include 'COMMON.INTERACT'
1057       include 'COMMON.IOUNITS'
1058       include 'COMMON.NAMES'
1059       integer i,j,inres
1060       do j=1,3
1061         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1062       enddo
1063       do i=nnt,nct-1
1064         do j=1,3
1065           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1066         enddo
1067       enddo
1068       do i=nnt,nct
1069         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1070           inres=i+nres
1071           do j=1,3
1072             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1073           enddo
1074         endif
1075       enddo
1076       return
1077       end
1078 c-----------------------------------------------------------------
1079       subroutine verlet1
1080 c Applying velocity Verlet algorithm - step 1 to coordinates
1081       implicit none
1082       include 'DIMENSIONS'
1083       include 'COMMON.CONTROL'
1084       include 'COMMON.VAR'
1085       include 'COMMON.MD'
1086 #ifdef FIVEDIAG
1087        include 'COMMON.LAGRANGE.5diag'
1088 #else
1089        include 'COMMON.LAGRANGE'
1090 #endif
1091       include 'COMMON.CHAIN'
1092       include 'COMMON.DERIV'
1093       include 'COMMON.GEO'
1094       include 'COMMON.LOCAL'
1095       include 'COMMON.INTERACT'
1096       include 'COMMON.IOUNITS'
1097       include 'COMMON.NAMES'
1098       double precision adt,adt2
1099       integer i,j,inres
1100 #ifdef DEBUG
1101       write (iout,*) "VELVERLET1 START: DC"
1102       do i=0,nres
1103         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1104      &   (dc(j,i+nres),j=1,3)
1105       enddo
1106 #endif
1107       do j=1,3
1108         adt=d_a_old(j,0)*d_time
1109         adt2=0.5d0*adt
1110         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1111         d_t_new(j,0)=d_t_old(j,0)+adt2
1112         d_t(j,0)=d_t_old(j,0)+adt
1113       enddo
1114       do i=nnt,nct-1
1115 C      SPYTAC ADAMA
1116 C       do i=0,nres
1117 #ifdef DEBUG
1118         write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
1119 #endif
1120         do j=1,3
1121           adt=d_a_old(j,i)*d_time
1122           adt2=0.5d0*adt
1123           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1124           d_t_new(j,i)=d_t_old(j,i)+adt2
1125           d_t(j,i)=d_t_old(j,i)+adt
1126         enddo
1127       enddo
1128       do i=nnt,nct
1129 C        do i=0,nres
1130         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1131           inres=i+nres
1132           do j=1,3
1133             adt=d_a_old(j,inres)*d_time
1134             adt2=0.5d0*adt
1135             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1136             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1137             d_t(j,inres)=d_t_old(j,inres)+adt
1138           enddo
1139         endif
1140       enddo
1141 #ifdef DEBUG
1142       write (iout,*) "VELVERLET1 END: DC"
1143       do i=0,nres
1144         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1145      &   (dc(j,i+nres),j=1,3)
1146       enddo
1147 #endif
1148       return
1149       end
1150 c---------------------------------------------------------------------
1151       subroutine verlet2
1152 c  Step 2 of the velocity Verlet algorithm: update velocities
1153       implicit none
1154       include 'DIMENSIONS'
1155       include 'COMMON.CONTROL'
1156       include 'COMMON.VAR'
1157       include 'COMMON.MD'
1158 #ifdef FIVEDIAG
1159        include 'COMMON.LAGRANGE.5diag'
1160 #else
1161        include 'COMMON.LAGRANGE'
1162 #endif
1163       include 'COMMON.CHAIN'
1164       include 'COMMON.DERIV'
1165       include 'COMMON.GEO'
1166       include 'COMMON.LOCAL'
1167       include 'COMMON.INTERACT'
1168       include 'COMMON.IOUNITS'
1169       include 'COMMON.NAMES'
1170       integer i,j,inres
1171       do j=1,3
1172         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1173       enddo
1174       do i=nnt,nct-1
1175         do j=1,3
1176           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1177         enddo
1178       enddo
1179       do i=nnt,nct
1180         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1181           inres=i+nres
1182           do j=1,3
1183             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1184           enddo
1185         endif
1186       enddo
1187       return
1188       end
1189 c-----------------------------------------------------------------
1190       subroutine sddir_precalc
1191 c Applying velocity Verlet algorithm - step 1 to coordinates
1192       implicit none
1193       include 'DIMENSIONS'
1194 #ifdef MPI
1195       include 'mpif.h'
1196 #endif
1197       include 'COMMON.CONTROL'
1198       include 'COMMON.VAR'
1199       include 'COMMON.MD'
1200 #ifdef FIVEDIAG
1201        include 'COMMON.LAGRANGE.5diag'
1202 #else
1203        include 'COMMON.LAGRANGE'
1204 #endif
1205 #ifndef LANG0
1206       include 'COMMON.LANGEVIN'
1207 #else
1208 #ifdef FIVEDIAG
1209       include 'COMMON.LANGEVIN.lang0.5diag'
1210 #else
1211       include 'COMMON.LANGEVIN.lang0'
1212 #endif
1213 #endif
1214       include 'COMMON.CHAIN'
1215       include 'COMMON.DERIV'
1216       include 'COMMON.GEO'
1217       include 'COMMON.LOCAL'
1218       include 'COMMON.INTERACT'
1219       include 'COMMON.IOUNITS'
1220       include 'COMMON.NAMES'
1221       include 'COMMON.TIME1'
1222       double precision time00
1223       double precision stochforcvec(MAXRES6)
1224       common /stochcalc/ stochforcvec
1225       integer i
1226 c
1227 c Compute friction and stochastic forces
1228 c
1229 #ifdef MPI
1230       time00=MPI_Wtime()
1231 #else
1232       time00=tcpu()
1233 #endif
1234       call friction_force
1235 c      write (iout,*) "After friction_force"
1236 #ifdef MPI
1237       time_fric=time_fric+MPI_Wtime()-time00
1238       time00=MPI_Wtime()
1239 #else
1240       time_fric=time_fric+tcpu()-time00
1241       time00=tcpu()
1242 #endif
1243       call stochastic_force(stochforcvec)
1244 c      write (iout,*) "After stochastic_force"
1245 #ifdef MPI
1246       time_stoch=time_stoch+MPI_Wtime()-time00
1247 #else
1248       time_stoch=time_stoch+tcpu()-time00
1249 #endif
1250 c
1251 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1252 c forces (d_as_work)
1253 c
1254 #ifdef FIVEDIAG
1255 c      write (iout,*) "friction accelerations"
1256       call fivediaginv_mult(dimen,fric_work, d_af_work)
1257 c      write (iout,*) "stochastic acceleratios"
1258       call fivediaginv_mult(dimen,stochforcvec, d_as_work)
1259 c      write (iout,*) "Leaving sddir_precalc"
1260 #else
1261       call ginv_mult(fric_work, d_af_work)
1262       call ginv_mult(stochforcvec, d_as_work)
1263 #endif
1264 #ifdef DEBUG
1265       write (iout,*) "d_af_work"
1266       write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
1267       write (iout,*) "d_as_work"
1268       write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
1269 #endif
1270       return
1271       end
1272 c---------------------------------------------------------------------
1273       subroutine sddir_verlet1
1274 c Applying velocity Verlet algorithm - step 1 to velocities
1275       implicit none
1276       include 'DIMENSIONS'
1277       include 'COMMON.CONTROL'
1278       include 'COMMON.VAR'
1279       include 'COMMON.MD'
1280 #ifdef FIVEDIAG
1281        include 'COMMON.LAGRANGE.5diag'
1282 #else
1283        include 'COMMON.LAGRANGE'
1284 #endif
1285 #ifndef LANG0
1286       include 'COMMON.LANGEVIN'
1287 #else
1288 #ifdef FIVEDIAG
1289       include 'COMMON.LANGEVIN.lang0.5diag'
1290 #else
1291       include 'COMMON.LANGEVIN.lang0'
1292 #endif
1293 #endif
1294       include 'COMMON.CHAIN'
1295       include 'COMMON.DERIV'
1296       include 'COMMON.GEO'
1297       include 'COMMON.LOCAL'
1298       include 'COMMON.INTERACT'
1299       include 'COMMON.IOUNITS'
1300       include 'COMMON.NAMES'
1301 c Revised 3/31/05 AL: correlation between random contributions to
1302 c position and velocity increments included.
1303       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1304       double precision adt,adt2
1305       integer i,j,ind,inres
1306 c
1307 c Add the contribution from BOTH friction and stochastic force to the
1308 c coordinates, but ONLY the contribution from the friction forces to velocities
1309 c
1310       do j=1,3
1311         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1312         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1313         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1314         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1315         d_t(j,0)=d_t_old(j,0)+adt
1316       enddo
1317       ind=3
1318       do i=nnt,nct-1
1319         do j=1,3
1320           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1321           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1322           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1323           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1324           d_t(j,i)=d_t_old(j,i)+adt
1325         enddo
1326         ind=ind+3
1327       enddo
1328       do i=nnt,nct
1329         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1330           inres=i+nres
1331           do j=1,3
1332             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1333             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1334             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1335             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1336             d_t(j,inres)=d_t_old(j,inres)+adt
1337           enddo
1338           ind=ind+3
1339         endif
1340       enddo
1341       return
1342       end
1343 c---------------------------------------------------------------------
1344       subroutine sddir_verlet2
1345 c  Calculating the adjusted velocities for accelerations
1346       implicit none
1347       include 'DIMENSIONS'
1348       include 'COMMON.CONTROL'
1349       include 'COMMON.VAR'
1350       include 'COMMON.MD'
1351 #ifdef FIVEDIAG
1352        include 'COMMON.LAGRANGE.5diag'
1353 #else
1354        include 'COMMON.LAGRANGE'
1355 #endif
1356 #ifndef LANG0
1357       include 'COMMON.LANGEVIN'
1358 #else
1359 #ifdef FIVEDIAG
1360       include 'COMMON.LANGEVIN.lang0.5diag'
1361 #else
1362       include 'COMMON.LANGEVIN.lang0'
1363 #endif
1364 #endif
1365       include 'COMMON.CHAIN'
1366       include 'COMMON.DERIV'
1367       include 'COMMON.GEO'
1368       include 'COMMON.LOCAL'
1369       include 'COMMON.INTERACT'
1370       include 'COMMON.IOUNITS'
1371       include 'COMMON.NAMES'
1372       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1373       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1374       integer i,j,inres,ind
1375 c Revised 3/31/05 AL: correlation between random contributions to
1376 c position and velocity increments included.
1377 c The correlation coefficients are calculated at low-friction limit.
1378 c Also, friction forces are now not calculated with new velocities.
1379
1380 c      call friction_force
1381       call stochastic_force(stochforcvec)
1382 c
1383 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1384 c forces (d_as_work)
1385 c
1386 #ifdef FIVEDIAG
1387       call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
1388 #else
1389       call ginv_mult(stochforcvec, d_as_work1)
1390 #endif
1391 c
1392 c Update velocities
1393 c
1394       do j=1,3
1395         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1396      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1397       enddo
1398       ind=3
1399       do i=nnt,nct-1
1400         do j=1,3
1401           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1402      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1403         enddo
1404         ind=ind+3
1405       enddo
1406       do i=nnt,nct
1407         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1408           inres=i+nres
1409           do j=1,3
1410             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1411      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1412      &       +cos60*d_as_work1(ind+j))*d_time
1413           enddo
1414           ind=ind+3
1415         endif
1416       enddo
1417       return
1418       end
1419 c---------------------------------------------------------------------
1420       subroutine max_accel
1421 c
1422 c Find the maximum difference in the accelerations of the the sites
1423 c at the beginning and the end of the time step.
1424 c
1425       implicit none
1426       include 'DIMENSIONS'
1427       include 'COMMON.CONTROL'
1428       include 'COMMON.VAR'
1429       include 'COMMON.MD'
1430 #ifdef FIVEDIAG
1431        include 'COMMON.LAGRANGE.5diag'
1432 #else
1433        include 'COMMON.LAGRANGE'
1434 #endif
1435       include 'COMMON.CHAIN'
1436       include 'COMMON.DERIV'
1437       include 'COMMON.GEO'
1438       include 'COMMON.LOCAL'
1439       include 'COMMON.INTERACT'
1440       include 'COMMON.IOUNITS'
1441       integer i,j
1442       double precision aux(3),accel(3),accel_old(3),dacc
1443       do j=1,3
1444 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1445          accel_old(j)=d_a_old(j,0)
1446          accel(j)=d_a(j,0)
1447       enddo
1448       amax=0.0d0
1449       do i=nnt,nct
1450 c Backbone
1451         if (i.lt.nct) then
1452 c 7/3/08 changed to asymmetric difference
1453           do j=1,3
1454 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1455             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1456             accel(j)=accel(j)+0.5d0*d_a(j,i)
1457 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1458             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1459               dacc=dabs(accel(j)-accel_old(j))
1460 c              write (iout,*) i,dacc
1461               if (dacc.gt.amax) amax=dacc
1462             endif
1463           enddo
1464         endif
1465       enddo
1466 c Side chains
1467       do j=1,3
1468 c        accel(j)=aux(j)
1469         accel_old(j)=d_a_old(j,0)
1470         accel(j)=d_a(j,0)
1471       enddo
1472       if (nnt.eq.2) then
1473         do j=1,3
1474           accel_old(j)=accel_old(j)+d_a_old(j,1)
1475           accel(j)=accel(j)+d_a(j,1)
1476         enddo
1477       endif
1478       do i=nnt,nct
1479         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1480           do j=1,3
1481 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1482             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1483             accel(j)=accel(j)+d_a(j,i+nres)
1484           enddo
1485         endif
1486         do j=1,3
1487 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1488           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1489             dacc=dabs(accel(j)-accel_old(j))
1490 c            write (iout,*) "side-chain",i,dacc
1491             if (dacc.gt.amax) amax=dacc
1492           endif
1493         enddo
1494         do j=1,3
1495           accel_old(j)=accel_old(j)+d_a_old(j,i)
1496           accel(j)=accel(j)+d_a(j,i)
1497 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1498         enddo
1499       enddo
1500       return
1501       end
1502 c---------------------------------------------------------------------
1503       subroutine predict_edrift(epdrift)
1504 c
1505 c Predict the drift of the potential energy
1506 c
1507       implicit none
1508       include 'DIMENSIONS'
1509       include 'COMMON.CONTROL'
1510       include 'COMMON.VAR'
1511       include 'COMMON.MD'
1512 #ifdef FIVEDIAG
1513        include 'COMMON.LAGRANGE.5diag'
1514 #else
1515        include 'COMMON.LAGRANGE'
1516 #endif
1517       include 'COMMON.CHAIN'
1518       include 'COMMON.DERIV'
1519       include 'COMMON.GEO'
1520       include 'COMMON.LOCAL'
1521       include 'COMMON.INTERACT'
1522       include 'COMMON.IOUNITS'
1523       include 'COMMON.MUCA'
1524       double precision epdrift,epdriftij
1525       integer i,j
1526 c Drift of the potential energy
1527       epdrift=0.0d0
1528       do i=nnt,nct
1529 c Backbone
1530         if (i.lt.nct) then
1531           do j=1,3
1532             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1533             if (lmuca) epdriftij=epdriftij*factor
1534 c            write (iout,*) "back",i,j,epdriftij
1535             if (epdriftij.gt.epdrift) epdrift=epdriftij
1536           enddo
1537         endif
1538 c Side chains
1539         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1540           do j=1,3
1541             epdriftij=
1542      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1543             if (lmuca) epdriftij=epdriftij*factor
1544 c            write (iout,*) "side",i,j,epdriftij
1545             if (epdriftij.gt.epdrift) epdrift=epdriftij
1546           enddo
1547         endif
1548       enddo
1549       epdrift=0.5d0*epdrift*d_time*d_time
1550 c      write (iout,*) "epdrift",epdrift
1551       return
1552       end
1553 c-----------------------------------------------------------------------
1554       subroutine verlet_bath
1555 c
1556 c  Coupling to the thermostat by using the Berendsen algorithm
1557 c
1558       implicit none
1559       include 'DIMENSIONS'
1560       include 'COMMON.CONTROL'
1561       include 'COMMON.VAR'
1562       include 'COMMON.MD'
1563 #ifdef FIVEDIAG
1564        include 'COMMON.LAGRANGE.5diag'
1565 #else
1566        include 'COMMON.LAGRANGE'
1567 #endif
1568 #ifdef LANG0
1569 #ifdef FIVEDIAG
1570       include 'COMMON.LANGEVIN.lang0.5diag'
1571 #else
1572       include 'COMMON.LANGEVIN.lang0'
1573 #endif
1574 #else
1575       include 'COMMON.LANGEVIN'
1576 #endif
1577       include 'COMMON.CHAIN'
1578       include 'COMMON.DERIV'
1579       include 'COMMON.GEO'
1580       include 'COMMON.LOCAL'
1581       include 'COMMON.INTERACT'
1582       include 'COMMON.IOUNITS'
1583       include 'COMMON.NAMES'
1584       double precision T_half,fact
1585       integer i,j ,inres
1586       T_half=2.0d0/(dimen3*Rb)*EK
1587       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1588 c      write(iout,*) "T_half", T_half
1589 c      write(iout,*) "EK", EK
1590 c      write(iout,*) "fact", fact
1591       do j=1,3
1592         d_t(j,0)=fact*d_t(j,0)
1593       enddo
1594       do i=nnt,nct-1
1595         do j=1,3
1596           d_t(j,i)=fact*d_t(j,i)
1597         enddo
1598       enddo
1599       do i=nnt,nct
1600         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1601           inres=i+nres
1602           do j=1,3
1603             d_t(j,inres)=fact*d_t(j,inres)
1604           enddo
1605         endif
1606       enddo
1607       return
1608       end
1609 c---------------------------------------------------------
1610       subroutine init_MD
1611 c  Set up the initial conditions of a MD simulation
1612       implicit none
1613       include 'DIMENSIONS'
1614 #ifdef MP
1615       include 'mpif.h'
1616       character*16 form
1617       integer IERROR,ERRCODE,error_msg,ierr,ierrcode
1618 #endif
1619       include 'COMMON.SETUP'
1620       include 'COMMON.CONTROL'
1621       include 'COMMON.VAR'
1622       include 'COMMON.MD'
1623 #ifdef FIVEDIAG
1624        include 'COMMON.LAGRANGE.5diag'
1625 #else
1626        include 'COMMON.LAGRANGE'
1627 #endif
1628       include 'COMMON.QRESTR'
1629 #ifndef LANG0
1630       include 'COMMON.LANGEVIN'
1631 #else
1632 #ifdef FIVEDIAG
1633       include 'COMMON.LANGEVIN.lang0.5diag'
1634 #else
1635 #ifdef FIVEDIAG
1636       include 'COMMON.LANGEVIN.lang0.5diag'
1637 #else
1638       include 'COMMON.LANGEVIN.lang0'
1639 #endif
1640 #endif
1641 #endif
1642       include 'COMMON.CHAIN'
1643       include 'COMMON.DERIV'
1644       include 'COMMON.GEO'
1645       include 'COMMON.LOCAL'
1646       include 'COMMON.INTERACT'
1647       include 'COMMON.IOUNITS'
1648       include 'COMMON.NAMES'
1649       include 'COMMON.REMD'
1650       include 'COMMON.TIME1'
1651 #ifdef LBFGS
1652       character*9 status
1653       integer niter
1654       common /lbfgstat/ status,niter,nfun
1655 #endif
1656       integer n_model_try,list_model_try(max_template),k
1657       double precision tt0
1658       real*8 energia_long(0:n_ene),
1659      &  energia_short(0:n_ene),vcm(3),incr(3)
1660       double precision cm(3),L(3),xv,sigv,lowb,highb
1661       double precision varia(maxvar),energia(0:n_ene)
1662       character*256 qstr
1663       integer ilen
1664       external ilen
1665       character*50 tytul
1666       logical file_exist
1667       common /gucio/ cm
1668       integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
1669      & i_model,itime
1670       integer iran_num
1671       double precision etot
1672       logical fail
1673       write (iout,*) "init_MD INDPDB",indpdb
1674       d_time0=d_time
1675 c      write(iout,*) "d_time", d_time
1676 c Compute the standard deviations of stochastic forces for Langevin dynamics
1677 c if the friction coefficients do not depend on surface area
1678       if (lang.gt.0 .and. .not.surfarea) then
1679         do i=nnt,nct-1
1680           stdforcp(i)=stdfp*dsqrt(gamp)
1681         enddo
1682         do i=nnt,nct
1683           if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1684      &                *dsqrt(gamsc(iabs(itype(i))))
1685         enddo
1686       endif
1687 c Open the pdb file for snapshotshots
1688 #ifdef MPI
1689       if(mdpdb) then
1690         if (ilen(tmpdir).gt.0)
1691      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1692      &      liczba(:ilen(liczba))//".pdb")
1693         open(ipdb,
1694      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1695      &  //".pdb")
1696       else
1697 #ifdef NOXDR
1698         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1699      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1700      &      liczba(:ilen(liczba))//".x")
1701         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1702      &  //".x"
1703 #else
1704         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1705      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1706      &      liczba(:ilen(liczba))//".cx")
1707         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1708      &  //".cx"
1709 #endif
1710       endif
1711 #else
1712       if(mdpdb) then
1713          if (ilen(tmpdir).gt.0)
1714      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1715          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1716       else
1717          if (ilen(tmpdir).gt.0)
1718      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1719          cartname=prefix(:ilen(prefix))//"_MD.cx"
1720       endif
1721 #endif
1722       if (usampl) then
1723         write (qstr,'(256(1h ))')
1724         ipos=1
1725         do i=1,nfrag
1726           iq = qinfrag(i,iset)*10
1727           iw = wfrag(i,iset)/100
1728           if (iw.gt.0) then
1729             if(me.eq.king.or..not.out1file)
1730      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1731             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1732             ipos=ipos+7
1733           endif
1734         enddo
1735         do i=1,npair
1736           iq = qinpair(i,iset)*10
1737           iw = wpair(i,iset)/100
1738           if (iw.gt.0) then
1739             if(me.eq.king.or..not.out1file)
1740      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1741             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1742             ipos=ipos+7
1743           endif
1744         enddo
1745 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1746 #ifdef NOXDR
1747 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1748 #else
1749 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1750 #endif
1751 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1752       endif
1753       icg=1
1754       write (iout,*) "REST ",rest
1755       if (rest) then
1756        if (restart1file) then
1757          if (me.eq.king)
1758      &     inquire(file=mremd_rst_name,exist=file_exist)
1759 #ifdef MPI
1760            write (*,*) me," Before broadcast: file_exist",file_exist
1761          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1762      &          IERR)
1763          write (*,*) me," After broadcast: file_exist",file_exist
1764 c        inquire(file=mremd_rst_name,exist=file_exist)
1765 #endif
1766         if(me.eq.king.or..not.out1file)
1767      &   write(iout,*) "Initial state read by master and distributed"
1768        else
1769          if (ilen(tmpdir).gt.0)
1770      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1771      &      //liczba(:ilen(liczba))//'.rst')
1772         inquire(file=rest2name,exist=file_exist)
1773        endif
1774        if(file_exist) then
1775          if(.not.restart1file) then
1776            if(me.eq.king.or..not.out1file)
1777      &      write(iout,*) "Initial state will be read from file ",
1778      &      rest2name(:ilen(rest2name))
1779            call readrst
1780          endif
1781          call rescale_weights(t_bath)
1782        else
1783         rest=.false.
1784         if(me.eq.king.or..not.out1file)then
1785          if (restart1file) then
1786           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1787      &       " does not exist"
1788          else
1789           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1790      &       " does not exist"
1791          endif
1792          write(iout,*) "Initial velocities randomly generated"
1793         endif
1794         call random_vel
1795         totT=0.0d0
1796         totTafm=totT
1797        endif
1798       else
1799 c Generate initial velocities
1800         if(me.eq.king.or..not.out1file)
1801      &   write(iout,*) "Initial velocities randomly generated"
1802         call random_vel
1803         totT=0.0d0
1804 CtotTafm is the variable for AFM time which eclipsed during
1805         totTafm=totT
1806       endif
1807 c      rest2name = prefix(:ilen(prefix))//'.rst'
1808       if(me.eq.king.or..not.out1file)then
1809        write (iout,*) "Initial velocities"
1810        do i=0,nres
1811          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1812      &   (d_t(j,i+nres),j=1,3)
1813        enddo
1814 c  Zeroing the total angular momentum of the system
1815        write(iout,*) "Calling the zero-angular
1816      & momentum subroutine"
1817       endif
1818       call inertia_tensor
1819 c  Getting the potential energy and forces and velocities and accelerations
1820       call vcm_vel(vcm)
1821       write (iout,*) "velocity of the center of the mass:"
1822       write (iout,*) (vcm(j),j=1,3)
1823       call flush(iout)
1824       do j=1,3
1825         d_t(j,0)=d_t(j,0)-vcm(j)
1826       enddo
1827 c Removing the velocity of the center of mass
1828       call vcm_vel(vcm)
1829       if(me.eq.king.or..not.out1file)then
1830        write (iout,*) "vcm right after adjustment:"
1831        write (iout,*) (vcm(j),j=1,3)
1832       endif
1833       call flush(iout)
1834       write (iout,*) "init_MD before initial structure REST ",rest
1835       if (.not.rest) then
1836   122   continue
1837         if (iranconf.ne.0) then
1838 c 8/22/17 AL Loop to produce a low-energy random conformation
1839           do iranmin=1,10
1840             call chainbuild
1841             if (overlapsc) then
1842               print *, 'Calling OVERLAP_SC'
1843               call overlap_sc(fail)
1844             endif
1845
1846             if (searchsc) then
1847               call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1848               print *,'SC_move',nft_sc,etot
1849               if(me.eq.king.or..not.out1file)
1850      &          write(iout,*) 'SC_move',nft_sc,etot
1851             endif
1852
1853             if (dccart) then
1854               if (me.eq.king.or..not.out1file) write(iout,*)
1855      &         'Minimizing random structure: Calling MINIM_DC'
1856                 call minim_dc(etot,iretcode,nfun)
1857             else
1858               call geom_to_var(nvar,varia)
1859               if (me.eq.king.or..not.out1file) write (iout,*)
1860      &          'Minimizing random structure: Calling MINIMIZE.'
1861               call minimize(etot,varia,iretcode,nfun)
1862               call var_to_geom(nvar,varia)
1863             endif
1864             if (me.eq.king.or..not.out1file)
1865      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1866             if (isnan(etot) .or. etot.gt.1.0d4) then
1867                write (iout,*) "Energy too large",etot,
1868      &          " trying another random conformation"
1869               do itrial=1,100
1870                 itmp=1
1871                 call gen_rand_conf(itmp,*30)
1872                 goto 40
1873    30           write (iout,*) 'Failed to generate random conformation',
1874      &            ', itrial=',itrial
1875                 write (*,*) 'Processor:',me,
1876      &            ' Failed to generate random conformation',
1877      &            ' itrial=',itrial
1878                 call intout
1879
1880 #ifdef AIX
1881                 call flush_(iout)
1882 #else
1883                 call flush(iout)
1884 #endif
1885               enddo
1886               write (iout,'(a,i3,a)') 'Processor:',me,
1887      &        ' error in generating random conformation.'
1888               write (*,'(a,i3,a)') 'Processor:',me,
1889      &        ' error in generating random conformation.'
1890               call flush(iout)
1891 #ifdef MPI
1892               call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1893 #else
1894               stop
1895 #endif
1896    40         continue
1897             else
1898               goto 44
1899             endif
1900           enddo
1901           write (iout,'(a,i3,a)') 'Processor:',me,
1902      &        ' failed to generate a low-energy random conformation.'
1903           write (*,'(a,i3,a)') 'Processor:',me,
1904      &        ' failed to generate a low-energy random conformation.'
1905             call flush(iout)
1906 #ifdef MPI
1907             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1908 #else
1909             stop
1910 #endif
1911    44     continue
1912         else if (preminim) then
1913           if (start_from_model) then
1914             n_model_try=0
1915             do while (fail .and. n_model_try.lt.constr_homology)
1916               do
1917                 i_model=iran_num(1,constr_homology)
1918                 do k=1,n_model_try
1919                   if (i_model.eq.list_model_try(k)) exit
1920                 enddo
1921                 if (k.gt.n_model_try) exit
1922               enddo
1923               n_model_try=n_model_try+1
1924               list_model_try(n_model_try)=i_model
1925               write (iout,*) 'starting from model ',i_model
1926               do i=1,2*nres
1927                 do j=1,3
1928                   c(j,i)=chomo(j,i,i_model)
1929                 enddo
1930               enddo
1931               call int_from_cart(.true.,.false.)
1932               call sc_loc_geom(.false.)
1933               do i=1,nres-1
1934                 do j=1,3
1935                   dc(j,i)=c(j,i+1)-c(j,i)
1936                   dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1937                 enddo
1938               enddo
1939               do i=2,nres-1
1940                 do j=1,3
1941                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1942                   dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1943                 enddo
1944               enddo
1945               if (me.eq.king.or..not.out1file) then
1946               write (iout,*) "Energies before removing overlaps"
1947               call etotal(energia(0))
1948               call enerprint(energia(0))
1949               endif
1950 ! Remove SC overlaps if requested
1951               if (overlapsc) then
1952                 write (iout,*) 'Calling OVERLAP_SC'
1953                 call overlap_sc(fail)
1954                 if (fail) then
1955                   write (iout,*)
1956      &            "Failed to remove overlap from model",i_model
1957                   cycle
1958                 endif
1959               endif
1960 #ifdef SEARCHSC
1961               if (me.eq.king.or..not.out1file) then
1962               write (iout,*) "Energies after removing overlaps"
1963               call etotal(energia(0))
1964               call enerprint(energia(0))
1965               endif
1966 ! Search for better SC rotamers if requested
1967               if (searchsc) then
1968                 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1969                 print *,'SC_move',nft_sc,etot
1970                 if (me.eq.king.or..not.out1file)
1971      &            write(iout,*) 'SC_move',nft_sc,etot
1972               endif
1973               call etotal(energia(0))
1974 #endif
1975             enddo
1976             if (n_model_try.gt.constr_homology) then
1977               write (iout,*)
1978      &    "All models have irreparable overlaps. Trying randoms starts."
1979               iranconf=1
1980               goto 122
1981             endif
1982           endif
1983 C 8/22/17 AL Minimize initial structure
1984           if (dccart) then
1985             if (me.eq.king.or..not.out1file) write(iout,*)
1986      &        'Minimizing initial PDB structure: Calling MINIM_DC'
1987             call minim_dc(etot,iretcode,nfun)
1988           else
1989             call geom_to_var(nvar,varia)
1990             if(me.eq.king.or..not.out1file) write (iout,*)
1991      &        'Minimizing initial PDB structure: Calling MINIMIZE.'
1992             call minimize(etot,varia,iretcode,nfun)
1993             call var_to_geom(nvar,varia)
1994 #ifdef LBFGS
1995             if (me.eq.king.or..not.out1file)
1996      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
1997             if(me.eq.king.or..not.out1file)
1998      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
1999 #else
2000             if (me.eq.king.or..not.out1file)
2001      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2002             if(me.eq.king.or..not.out1file)
2003      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2004 #endif
2005           endif
2006         endif
2007       endif ! .not. rest
2008       call chainbuild_cart
2009       call kinetic(EK)
2010       if (tbf) then
2011         call verlet_bath
2012       endif
2013       kinetic_T=2.0d0/(dimen3*Rb)*EK
2014       if(me.eq.king.or..not.out1file)then
2015        call cartprint
2016        call intout
2017       endif
2018 #ifdef MPI
2019       tt0=MPI_Wtime()
2020 #else
2021       tt0=tcpu()
2022 #endif
2023       call zerograd
2024       call etotal(potEcomp)
2025       if (large) call enerprint(potEcomp)
2026 #ifdef TIMING_ENE
2027 #ifdef MPI
2028       t_etotal=t_etotal+MPI_Wtime()-tt0
2029 #else
2030       t_etotal=t_etotal+tcpu()-tt0
2031 #endif
2032 #endif
2033       potE=potEcomp(0)
2034 c      write (iout,*) "PotE-homology",potE-potEcomp(27)
2035       call cartgrad
2036       call lagrangian
2037       call max_accel
2038       if (amax*d_time .gt. dvmax) then
2039         d_time=d_time*dvmax/amax
2040         if(me.eq.king.or..not.out1file) write (iout,*)
2041      &   "Time step reduced to",d_time,
2042      &   " because of too large initial acceleration."
2043       endif
2044       if(me.eq.king.or..not.out1file)then
2045        write(iout,*) "Potential energy and its components"
2046        call enerprint(potEcomp)
2047 c       write(iout,*) (potEcomp(i),i=0,n_ene)
2048       endif
2049       potE=potEcomp(0)-potEcomp(27)
2050 c      write (iout,*) "PotE-homology",potE
2051       totE=EK+potE
2052       itime=0
2053       if (ntwe.ne.0) call statout(itime)
2054       if(me.eq.king.or..not.out1file)
2055      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2056      &   " Kinetic energy",EK," Potential energy",potE,
2057      &   " Total energy",totE," Maximum acceleration ",
2058      &   amax
2059       if (large) then
2060         write (iout,*) "Initial coordinates"
2061         do i=1,nres
2062           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2063      &    (c(j,i+nres),j=1,3)
2064         enddo
2065         write (iout,*) "Initial dC"
2066         do i=0,nres
2067           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2068      &    (dc(j,i+nres),j=1,3)
2069         enddo
2070         write (iout,*) "Initial velocities"
2071         write (iout,"(13x,' backbone ',23x,' side chain')")
2072         do i=0,nres
2073           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2074      &    (d_t(j,i+nres),j=1,3)
2075         enddo
2076         write (iout,*) "Initial accelerations"
2077         do i=0,nres
2078 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2079           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2080      &    (d_a(j,i+nres),j=1,3)
2081         enddo
2082       endif
2083       do i=0,2*nres
2084         do j=1,3
2085           dc_old(j,i)=dc(j,i)
2086           d_t_old(j,i)=d_t(j,i)
2087           d_a_old(j,i)=d_a(j,i)
2088         enddo
2089 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2090       enddo
2091       if (RESPA) then
2092 #ifdef MPI
2093         tt0 =MPI_Wtime()
2094 #else
2095         tt0 = tcpu()
2096 #endif
2097         call zerograd
2098         call etotal_short(energia_short)
2099         if (large) call enerprint(potEcomp)
2100 #ifdef TIMING_ENE
2101 #ifdef MPI
2102         t_eshort=t_eshort+MPI_Wtime()-tt0
2103 #else
2104         t_eshort=t_eshort+tcpu()-tt0
2105 #endif
2106 #endif
2107         call cartgrad
2108         call lagrangian
2109         if(.not.out1file .and. large) then
2110           write (iout,*) "energia_long",energia_long(0),
2111      &     " energia_short",energia_short(0),
2112      &     " total",energia_long(0)+energia_short(0)
2113           write (iout,*) "Initial fast-force accelerations"
2114           do i=0,nres
2115             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2116      &      (d_a(j,i+nres),j=1,3)
2117           enddo
2118         endif
2119 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2120         do i=0,2*nres
2121           do j=1,3
2122             d_a_short(j,i)=d_a(j,i)
2123           enddo
2124         enddo
2125 #ifdef MPI
2126         tt0=MPI_Wtime()
2127 #else
2128         tt0=tcpu()
2129 #endif
2130         call zerograd
2131         call etotal_long(energia_long)
2132         if (large) call enerprint(potEcomp)
2133 #ifdef TIMING_ENE
2134 #ifdef MPI
2135         t_elong=t_elong+MPI_Wtime()-tt0
2136 #else
2137         t_elong=t_elong+tcpu()-tt0
2138 #endif
2139 #endif
2140         call cartgrad
2141         call lagrangian
2142         if(.not.out1file .and. large) then
2143           write (iout,*) "energia_long",energia_long(0)
2144           write (iout,*) "Initial slow-force accelerations"
2145           do i=0,nres
2146             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2147      &      (d_a(j,i+nres),j=1,3)
2148           enddo
2149         endif
2150 #ifdef MPI
2151         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2152 #else
2153         t_enegrad=t_enegrad+tcpu()-tt0
2154 #endif
2155       endif
2156       return
2157       end
2158 c-----------------------------------------------------------
2159       subroutine random_vel
2160       implicit none
2161       include 'DIMENSIONS'
2162       include 'COMMON.CONTROL'
2163       include 'COMMON.VAR'
2164       include 'COMMON.MD'
2165 #ifdef FIVEDIAG
2166        include 'COMMON.LAGRANGE.5diag'
2167 #else
2168        include 'COMMON.LAGRANGE'
2169 #endif
2170 #ifndef LANG0
2171       include 'COMMON.LANGEVIN'
2172 #else
2173 #ifdef FIVEDIAG
2174       include 'COMMON.LANGEVIN.lang0.5diag'
2175 #else
2176       include 'COMMON.LANGEVIN.lang0'
2177 #endif
2178 #endif
2179       include 'COMMON.CHAIN'
2180       include 'COMMON.DERIV'
2181       include 'COMMON.GEO'
2182       include 'COMMON.LOCAL'
2183       include 'COMMON.INTERACT'
2184       include 'COMMON.IOUNITS'
2185       include 'COMMON.NAMES'
2186       include 'COMMON.TIME1'
2187       double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2188       integer i,ii,j,k,l,ind
2189       double precision anorm_distr
2190       logical lprn /.true./
2191 #ifdef FIVEDIAG
2192       integer ichain,n,innt,inct,ibeg,ierr
2193       double precision work(8*maxres6)
2194       integer iwork(maxres6)
2195       double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2196      & Gvec(maxres2_chain,maxres2_chain)
2197       common /przechowalnia/Ghalf,Geigen,Gvec
2198 #ifdef DEBUG
2199       double precision inertia(maxres2_chain,maxres2_chain)
2200 #endif
2201 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2202 c First generate velocities in the eigenspace of the G matrix
2203 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2204 c      call flush(iout)
2205 #ifdef DEBUG
2206       write (iout,*) "Random_vel, fivediag"
2207 #endif
2208       d_t=0.0d0
2209       Ek2=0.0d0
2210       EK=0.0d0
2211       Ek3=0.0d0
2212       do ichain=1,nchain
2213         ind=0
2214         ghalf=0.0d0
2215         n=dimen_chain(ichain)
2216         innt=iposd_chain(ichain)
2217         inct=innt+n-1
2218 #ifdef DEBUG
2219         write (iout,*) "Chain",ichain," n",n," start",innt
2220         do i=innt,inct
2221           if (i.lt.inct-1) then
2222             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2223      &         DU2orig(i)
2224           else if (i.eq.inct-1) then
2225             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2226           else
2227             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2228           endif
2229         enddo
2230 #endif
2231         ghalf(ind+1)=dmorig(innt)
2232         ghalf(ind+2)=du1orig(innt)
2233         ghalf(ind+3)=dmorig(innt+1)
2234         ind=ind+3
2235         do i=3,n
2236           ind=ind+i-3
2237 c          write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2238 c     &       " indu1",innt+i-1," indm",innt+i
2239           ghalf(ind+1)=du2orig(innt-1+i-2)
2240           ghalf(ind+2)=du1orig(innt-1+i-1)
2241           ghalf(ind+3)=dmorig(innt-1+i)
2242 c          write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2243 c     &       "DU1",innt-1+i-1,"DM ",innt-1+i
2244           ind=ind+3
2245         enddo
2246 #ifdef DEBUG
2247         ind=0
2248         do i=1,n
2249           do j=1,i
2250             ind=ind+1
2251             inertia(i,j)=ghalf(ind)
2252             inertia(j,i)=ghalf(ind)
2253           enddo
2254         enddo
2255 #endif
2256 #ifdef DEBUG
2257         write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2258         write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2259         call matoutr(n,ghalf)
2260 #endif
2261         call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2262         if (large) then
2263           write (iout,'(//a,i3)')
2264      &    "Eigenvectors and eigenvalues of the G matrix chain",ichain
2265           call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2266         endif
2267 #ifdef DIAGCHECK
2268 c check diagonalization
2269         do i=1,n
2270           do j=1,n
2271             aux=0.0d0
2272             do k=1,n
2273               do l=1,n
2274                 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2275               enddo
2276             enddo
2277             if (i.eq.j) then
2278               write (iout,*) i,j,aux,geigen(i)
2279             else
2280               write (iout,*) i,j,aux
2281             endif
2282           enddo
2283         enddo
2284 #endif
2285         xv=0.0d0
2286         ii=0
2287         do i=1,n
2288           do k=1,3
2289             ii=ii+1
2290             sigv=dsqrt((Rb*t_bath)/geigen(i))
2291             lowb=-5*sigv
2292             highb=5*sigv
2293             d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2294             EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2295 c            write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2296 c     &      " d_t_work_new",d_t_work_new(ii)
2297           enddo
2298         enddo
2299         do k=1,3
2300           do i=1,n
2301             ind=(i-1)*3+k
2302             d_t_work(ind)=0.0d0
2303             do j=1,n
2304               d_t_work(ind)=d_t_work(ind)
2305      &                +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2306             enddo
2307 c            write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2308 c            call flush(iout)
2309           enddo
2310         enddo
2311 #ifdef DEBUG
2312         aux=0.0d0
2313         do k=1,3
2314           do i=1,n
2315             do j=1,n
2316             aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2317             enddo
2318           enddo
2319         enddo
2320         Ek3=Ek3+aux/2
2321 #endif
2322 c Transfer to the d_t vector
2323         innt=chain_border(1,ichain)
2324         inct=chain_border(2,ichain)
2325         ind=0
2326 c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
2327         do i=innt,inct
2328           do j=1,3
2329             ind=ind+1
2330             d_t(j,i)=d_t_work(ind)
2331           enddo
2332           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2333             do j=1,3
2334               ind=ind+1
2335               d_t(j,i+nres)=d_t_work(ind)
2336             enddo
2337           endif
2338         enddo
2339       enddo
2340       if (large) then
2341         write (iout,*)
2342         write (iout,*) "Random velocities in the Calpha,SC space"
2343         do i=1,nres
2344           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2345      &    restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2346         enddo
2347       endif
2348       call kinetic_CASC(Ek1)
2349 !
2350 ! Transform the velocities to virtual-bond space
2351 !
2352 #define WLOS
2353 #ifdef WLOS
2354       do i=1,nres
2355         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2356           do j=1,3
2357             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2358           enddo
2359         else
2360           do j=1,3
2361             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2362             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2363           enddo
2364         end if
2365       enddo
2366       d_t(:,nct)=0.0d0
2367 c      d_a(:,0)=d_a(:,1)
2368 c      d_a(:,1)=0.0d0
2369 c      write (iout,*) "Shifting accelerations"
2370       do ichain=1,nchain
2371 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
2372 c     &     chain_border1(1,ichain)
2373         d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2374         d_t(:,chain_border1(1,ichain))=0.0d0
2375       enddo
2376 c      write (iout,*) "Adding accelerations"
2377       do ichain=2,nchain
2378 c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2379 c     &   chain_border(2,ichain-1)
2380         d_t(:,chain_border1(1,ichain)-1)=
2381      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2382         d_t(:,chain_border(2,ichain-1))=0.0d0
2383       enddo
2384       do ichain=2,nchain
2385         write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2386      &   chain_border(2,ichain-1)
2387         d_t(:,chain_border1(1,ichain)-1)=
2388      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2389         d_t(:,chain_border(2,ichain-1))=0.0d0
2390       enddo
2391 #else
2392       ibeg=0
2393 c      do j=1,3
2394 c        d_t(j,0)=d_t(j,nnt)
2395 c      enddo
2396       do ichain=1,nchain
2397       innt=chain_border(1,ichain)
2398       inct=chain_border(2,ichain)
2399 c      write (iout,*) "ichain",ichain," innt",innt," inct",inct
2400 c      write (iout,*) "ibeg",ibeg
2401       do j=1,3
2402         d_t(j,ibeg)=d_t(j,innt)
2403       enddo
2404       ibeg=inct+1
2405       do i=innt,inct
2406         if (iabs(itype(i).eq.10)) then
2407 c          write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2408           do j=1,3
2409             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2410           enddo
2411         else
2412           do j=1,3
2413             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2414             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2415           enddo
2416         end if
2417       enddo
2418       enddo
2419 #endif
2420       if (large) then
2421         write (iout,*)
2422         write (iout,*)
2423      &    "Random velocities in the virtual-bond-vector space"
2424         write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2425         do i=1,nres
2426           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2427      &     restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2428         enddo
2429         write (iout,*)
2430         write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2431      &   Ek
2432         write (iout,*)
2433      &   "Kinetic temperatures from inertia matrix eigenvalues",
2434      &   2*Ek/(3*dimen*Rb)
2435 #ifdef DEBUG
2436         write (iout,*) "Kinetic energy from inertia matrix",Ek3
2437         write (iout,*) "Kinetic temperatures from inertia",
2438      &   2*Ek3/(3*dimen*Rb)
2439 #endif
2440         write (iout,*) "Kinetic energy from velocities in CA-SC space",
2441      &   Ek1
2442         write (iout,*)
2443      &   "Kinetic temperatures from velovities in CA-SC space",
2444      &   2*Ek1/(3*dimen*Rb)
2445         call kinetic(Ek1)
2446         write (iout,*)
2447      &   "Kinetic energy from virtual-bond-vector velocities",Ek1
2448         write (iout,*)
2449      &   "Kinetic temperature from virtual-bond-vector velocities ",
2450      &   2*Ek1/(dimen3*Rb)
2451       endif
2452 #else
2453       xv=0.0d0
2454       ii=0
2455       do i=1,dimen
2456         do k=1,3
2457           ii=ii+1
2458           sigv=dsqrt((Rb*t_bath)/geigen(i))
2459           lowb=-5*sigv
2460           highb=5*sigv
2461           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2462
2463 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2464 c     &      " d_t_work_new",d_t_work_new(ii)
2465         enddo
2466       enddo
2467 C       if (SELFGUIDE.gt.0) then
2468 C       distance=0.0
2469 C       do j=1,3
2470 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2471 C       distance=distance+vec_afm(j)**2
2472 C       enddo
2473 C       distance=dsqrt(distance)
2474 C       do j=1,3
2475 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2476 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2477 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2478 C     &    d_t_work_new(j+(afmend-1)*3)
2479 C       enddo
2480
2481 C       endif
2482
2483 c diagnostics
2484 c      Ek1=0.0d0
2485 c      ii=0
2486 c      do i=1,dimen
2487 c        do k=1,3
2488 c          ii=ii+1
2489 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2490 c        enddo
2491 c      enddo
2492 c      write (iout,*) "Ek from eigenvectors",Ek1
2493 c end diagnostics
2494 c Transform velocities to UNRES coordinate space
2495       do k=0,2
2496         do i=1,dimen
2497           ind=(i-1)*3+k+1
2498           d_t_work(ind)=0.0d0
2499           do j=1,dimen
2500             d_t_work(ind)=d_t_work(ind)
2501      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2502           enddo
2503 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2504 c          call flush(iout)
2505         enddo
2506       enddo
2507 c Transfer to the d_t vector
2508       do j=1,3
2509         d_t(j,0)=d_t_work(j)
2510       enddo
2511       ind=3
2512       do i=nnt,nct-1
2513         do j=1,3
2514           ind=ind+1
2515           d_t(j,i)=d_t_work(ind)
2516         enddo
2517       enddo
2518       do i=nnt,nct
2519         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2520           do j=1,3
2521             ind=ind+1
2522             d_t(j,i+nres)=d_t_work(ind)
2523           enddo
2524         endif
2525       enddo
2526 c      call kinetic(EK)
2527 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2528 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2529 c      call flush(iout)
2530 #endif
2531       return
2532       end
2533 #ifndef LANG0
2534 c-----------------------------------------------------------
2535       subroutine sd_verlet_p_setup
2536 c Sets up the parameters of stochastic Verlet algorithm
2537       implicit none
2538       include 'DIMENSIONS'
2539 #ifdef MPI
2540       include 'mpif.h'
2541 #endif
2542       include 'COMMON.CONTROL'
2543       include 'COMMON.VAR'
2544       include 'COMMON.MD'
2545 #ifndef LANG0
2546       include 'COMMON.LANGEVIN'
2547 #else
2548 #ifdef FIVEDIAG
2549       include 'COMMON.LANGEVIN.lang0.5diag'
2550 #else
2551       include 'COMMON.LANGEVIN.lang0'
2552 #endif
2553 #endif
2554       include 'COMMON.CHAIN'
2555       include 'COMMON.DERIV'
2556       include 'COMMON.GEO'
2557       include 'COMMON.LOCAL'
2558       include 'COMMON.INTERACT'
2559       include 'COMMON.IOUNITS'
2560       include 'COMMON.NAMES'
2561       include 'COMMON.TIME1'
2562       double precision emgdt(MAXRES6),
2563      & pterm,vterm,rho,rhoc,vsig,
2564      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2565      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2566      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2567       logical lprn /.false./
2568       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2569       double precision ktm
2570 #ifdef MPI
2571       tt0 = MPI_Wtime()
2572 #else
2573       tt0 = tcpu()
2574 #endif
2575 c
2576 c AL 8/17/04 Code adapted from tinker
2577 c
2578 c Get the frictional and random terms for stochastic dynamics in the
2579 c eigenspace of mass-scaled UNRES friction matrix
2580 c
2581       do i = 1, dimen
2582             gdt = fricgam(i) * d_time
2583 c
2584 c Stochastic dynamics reduces to simple MD for zero friction
2585 c
2586             if (gdt .le. zero) then
2587                pfric_vec(i) = 1.0d0
2588                vfric_vec(i) = d_time
2589                afric_vec(i) = 0.5d0 * d_time * d_time
2590                prand_vec(i) = 0.0d0
2591                vrand_vec1(i) = 0.0d0
2592                vrand_vec2(i) = 0.0d0
2593 c
2594 c Analytical expressions when friction coefficient is large
2595 c
2596             else
2597                if (gdt .ge. gdt_radius) then
2598                   egdt = dexp(-gdt)
2599                   pfric_vec(i) = egdt
2600                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2601                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2602                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2603                   vterm = 1.0d0 - egdt**2
2604                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2605 c
2606 c Use series expansions when friction coefficient is small
2607 c
2608                else
2609                   gdt2 = gdt * gdt
2610                   gdt3 = gdt * gdt2
2611                   gdt4 = gdt2 * gdt2
2612                   gdt5 = gdt2 * gdt3
2613                   gdt6 = gdt3 * gdt3
2614                   gdt7 = gdt3 * gdt4
2615                   gdt8 = gdt4 * gdt4
2616                   gdt9 = gdt4 * gdt5
2617                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2618      &                          - gdt5/120.0d0 + gdt6/720.0d0
2619      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2620      &                          - gdt9/362880.0d0) / fricgam(i)**2
2621                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2622                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2623                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2624      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2625      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2626      &                       + 127.0d0*gdt9/90720.0d0
2627                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2628      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2629      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2630      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2631                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2632      &                       - 17.0d0*gdt2/1280.0d0
2633      &                       + 17.0d0*gdt3/6144.0d0
2634      &                       + 40967.0d0*gdt4/34406400.0d0
2635      &                       - 57203.0d0*gdt5/275251200.0d0
2636      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2637                end if
2638 c
2639 c Compute the scaling factors of random terms for the nonzero friction case
2640 c
2641                ktm = 0.5d0*d_time/fricgam(i)
2642                psig = dsqrt(ktm*pterm) / fricgam(i)
2643                vsig = dsqrt(ktm*vterm)
2644                rhoc = dsqrt(1.0d0 - rho*rho)
2645                prand_vec(i) = psig
2646                vrand_vec1(i) = vsig * rho
2647                vrand_vec2(i) = vsig * rhoc
2648             end if
2649       end do
2650       if (lprn) then
2651       write (iout,*)
2652      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2653      &  " vrand_vec2"
2654       do i=1,dimen
2655         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2656      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2657       enddo
2658       endif
2659 c
2660 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2661 c
2662       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2663       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2664       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2665       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2666       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2667       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2668 #ifdef MPI
2669       t_sdsetup=t_sdsetup+MPI_Wtime()
2670 #else
2671       t_sdsetup=t_sdsetup+tcpu()-tt0
2672 #endif
2673       return
2674       end
2675 c-------------------------------------------------------------
2676       subroutine eigtransf1(n,ndim,ab,d,c)
2677       implicit none
2678       integer n,ndim
2679       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2680       integer i,j,k
2681       do i=1,n
2682         do j=1,n
2683           c(i,j)=0.0d0
2684           do k=1,n
2685             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2686           enddo
2687         enddo
2688       enddo
2689       return
2690       end
2691 c-------------------------------------------------------------
2692       subroutine eigtransf(n,ndim,a,b,d,c)
2693       implicit none
2694       integer n,ndim
2695       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2696       integer i,j,k
2697       do i=1,n
2698         do j=1,n
2699           c(i,j)=0.0d0
2700           do k=1,n
2701             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2702           enddo
2703         enddo
2704       enddo
2705       return
2706       end
2707 c-------------------------------------------------------------
2708       subroutine sd_verlet1
2709 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2710       implicit none
2711       include 'DIMENSIONS'
2712       include 'COMMON.CONTROL'
2713       include 'COMMON.VAR'
2714       include 'COMMON.MD'
2715 #ifdef FIVEDIAG
2716        include 'COMMON.LAGRANGE.5diag'
2717 #else
2718        include 'COMMON.LAGRANGE'
2719 #endif
2720 #ifndef LANG0
2721       include 'COMMON.LANGEVIN'
2722 #else
2723 #ifdef FIVEDIAG
2724       include 'COMMON.LANGEVIN.lang0.5diag'
2725 #else
2726       include 'COMMON.LANGEVIN.lang0'
2727 #endif
2728 #endif
2729       include 'COMMON.CHAIN'
2730       include 'COMMON.DERIV'
2731       include 'COMMON.GEO'
2732       include 'COMMON.LOCAL'
2733       include 'COMMON.INTERACT'
2734       include 'COMMON.IOUNITS'
2735       include 'COMMON.NAMES'
2736       double precision stochforcvec(MAXRES6)
2737       common /stochcalc/ stochforcvec
2738       logical lprn /.false./
2739       integer i,j,ind,inres
2740
2741 c      write (iout,*) "dc_old"
2742 c      do i=0,nres
2743 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2744 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2745 c      enddo
2746       do j=1,3
2747         dc_work(j)=dc_old(j,0)
2748         d_t_work(j)=d_t_old(j,0)
2749         d_a_work(j)=d_a_old(j,0)
2750       enddo
2751       ind=3
2752       do i=nnt,nct-1
2753         do j=1,3
2754           dc_work(ind+j)=dc_old(j,i)
2755           d_t_work(ind+j)=d_t_old(j,i)
2756           d_a_work(ind+j)=d_a_old(j,i)
2757         enddo
2758         ind=ind+3
2759       enddo
2760       do i=nnt,nct
2761         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2762           do j=1,3
2763             dc_work(ind+j)=dc_old(j,i+nres)
2764             d_t_work(ind+j)=d_t_old(j,i+nres)
2765             d_a_work(ind+j)=d_a_old(j,i+nres)
2766           enddo
2767           ind=ind+3
2768         endif
2769       enddo
2770 #ifndef LANG0
2771       if (lprn) then
2772       write (iout,*)
2773      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2774      &  " vrand_mat2"
2775       do i=1,dimen
2776         do j=1,dimen
2777           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2778      &      vfric_mat(i,j),afric_mat(i,j),
2779      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2780         enddo
2781       enddo
2782       endif
2783       do i=1,dimen
2784         ddt1=0.0d0
2785         ddt2=0.0d0
2786         do j=1,dimen
2787           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2788      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2789           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2790           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2791         enddo
2792         d_t_work_new(i)=ddt1+0.5d0*ddt2
2793         d_t_work(i)=ddt1+ddt2
2794       enddo
2795 #endif
2796       do j=1,3
2797         dc(j,0)=dc_work(j)
2798         d_t(j,0)=d_t_work(j)
2799       enddo
2800       ind=3
2801       do i=nnt,nct-1
2802         do j=1,3
2803           dc(j,i)=dc_work(ind+j)
2804           d_t(j,i)=d_t_work(ind+j)
2805         enddo
2806         ind=ind+3
2807       enddo
2808       do i=nnt,nct
2809         if (itype(i).ne.10) then
2810           inres=i+nres
2811           do j=1,3
2812             dc(j,inres)=dc_work(ind+j)
2813             d_t(j,inres)=d_t_work(ind+j)
2814           enddo
2815           ind=ind+3
2816         endif
2817       enddo
2818       return
2819       end
2820 c--------------------------------------------------------------------------
2821       subroutine sd_verlet2
2822 c  Calculating the adjusted velocities for accelerations
2823       implicit none
2824       include 'DIMENSIONS'
2825       include 'COMMON.CONTROL'
2826       include 'COMMON.VAR'
2827       include 'COMMON.MD'
2828 #ifdef FIVEDIAG
2829        include 'COMMON.LAGRANGE.5diag'
2830 #else
2831        include 'COMMON.LAGRANGE'
2832 #endif
2833 #ifndef LANG0
2834       include 'COMMON.LANGEVIN'
2835 #else
2836 #ifdef FIVEDIAG
2837       include 'COMMON.LANGEVIN.lang0.5diag'
2838 #else
2839       include 'COMMON.LANGEVIN.lang0'
2840 #endif
2841 #endif
2842       include 'COMMON.CHAIN'
2843       include 'COMMON.DERIV'
2844       include 'COMMON.GEO'
2845       include 'COMMON.LOCAL'
2846       include 'COMMON.INTERACT'
2847       include 'COMMON.IOUNITS'
2848       include 'COMMON.NAMES'
2849       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2850       common /stochcalc/ stochforcvec
2851       integer i,j,ind,inres
2852 c
2853 c Compute the stochastic forces which contribute to velocity change
2854 c
2855       call stochastic_force(stochforcvecV)
2856
2857 #ifndef LANG0
2858       do i=1,dimen
2859         ddt1=0.0d0
2860         ddt2=0.0d0
2861         do j=1,dimen
2862           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2863           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2864      &     vrand_mat2(i,j)*stochforcvecV(j)
2865         enddo
2866         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2867       enddo
2868 #endif
2869       do j=1,3
2870         d_t(j,0)=d_t_work(j)
2871       enddo
2872       ind=3
2873       do i=nnt,nct-1
2874         do j=1,3
2875           d_t(j,i)=d_t_work(ind+j)
2876         enddo
2877         ind=ind+3
2878       enddo
2879       do i=nnt,nct
2880         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2881           inres=i+nres
2882           do j=1,3
2883             d_t(j,inres)=d_t_work(ind+j)
2884           enddo
2885           ind=ind+3
2886         endif
2887       enddo
2888       return
2889       end
2890 c-----------------------------------------------------------
2891       subroutine sd_verlet_ciccotti_setup
2892 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2893 c version
2894       implicit none
2895       include 'DIMENSIONS'
2896 #ifdef MPI
2897       include 'mpif.h'
2898 #endif
2899       include 'COMMON.CONTROL'
2900       include 'COMMON.VAR'
2901       include 'COMMON.MD'
2902 #ifdef FIVEDIAG
2903       include 'COMMON.LAGRANGE.5diag'
2904 #else
2905       include 'COMMON.LAGRANGE'
2906 #endif
2907       include 'COMMON.LANGEVIN'
2908       include 'COMMON.CHAIN'
2909       include 'COMMON.DERIV'
2910       include 'COMMON.GEO'
2911       include 'COMMON.LOCAL'
2912       include 'COMMON.INTERACT'
2913       include 'COMMON.IOUNITS'
2914       include 'COMMON.NAMES'
2915       include 'COMMON.TIME1'
2916       double precision emgdt(MAXRES6),
2917      & pterm,vterm,rho,rhoc,vsig,
2918      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2919      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2920      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2921       logical lprn /.false./
2922       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2923       double precision ktm
2924       integer i
2925 #ifdef MPI
2926       tt0 = MPI_Wtime()
2927 #else
2928       tt0 = tcpu()
2929 #endif
2930 c
2931 c AL 8/17/04 Code adapted from tinker
2932 c
2933 c Get the frictional and random terms for stochastic dynamics in the
2934 c eigenspace of mass-scaled UNRES friction matrix
2935 c
2936       do i = 1, dimen
2937             write (iout,*) "i",i," fricgam",fricgam(i)
2938             gdt = fricgam(i) * d_time
2939 c
2940 c Stochastic dynamics reduces to simple MD for zero friction
2941 c
2942             if (gdt .le. zero) then
2943                pfric_vec(i) = 1.0d0
2944                vfric_vec(i) = d_time
2945                afric_vec(i) = 0.5d0*d_time*d_time
2946                prand_vec(i) = afric_vec(i)
2947                vrand_vec2(i) = vfric_vec(i)
2948 c
2949 c Analytical expressions when friction coefficient is large
2950 c
2951             else
2952                egdt = dexp(-gdt)
2953                pfric_vec(i) = egdt
2954                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2955                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2956                prand_vec(i) = afric_vec(i)
2957                vrand_vec2(i) = vfric_vec(i)
2958 c
2959 c Compute the scaling factors of random terms for the nonzero friction case
2960 c
2961 c               ktm = 0.5d0*d_time/fricgam(i)
2962 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2963 c               vsig = dsqrt(ktm*vterm)
2964 c               prand_vec(i) = psig*afric_vec(i)
2965 c               vrand_vec2(i) = vsig*vfric_vec(i)
2966             end if
2967       end do
2968       if (lprn) then
2969       write (iout,*)
2970      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2971      &  " vrand_vec2"
2972       do i=1,dimen
2973         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2974      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2975       enddo
2976       endif
2977 c
2978 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2979 c
2980       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2981       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2982       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2983       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2984       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2985 #ifdef MPI
2986       t_sdsetup=t_sdsetup+MPI_Wtime()
2987 #else
2988       t_sdsetup=t_sdsetup+tcpu()-tt0
2989 #endif
2990       return
2991       end
2992 c-------------------------------------------------------------
2993       subroutine sd_verlet1_ciccotti
2994 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2995       implicit none
2996       include 'DIMENSIONS'
2997 #ifdef MPI
2998       include 'mpif.h'
2999 #endif
3000       include 'COMMON.CONTROL'
3001       include 'COMMON.VAR'
3002       include 'COMMON.MD'
3003 #ifdef FIVEDIAG
3004       include 'COMMON.LAGRANGE.5diag'
3005 #else
3006       include 'COMMON.LAGRANGE'
3007 #endif
3008       include 'COMMON.LANGEVIN'
3009       include 'COMMON.CHAIN'
3010       include 'COMMON.DERIV'
3011       include 'COMMON.GEO'
3012       include 'COMMON.LOCAL'
3013       include 'COMMON.INTERACT'
3014       include 'COMMON.IOUNITS'
3015       include 'COMMON.NAMES'
3016       double precision stochforcvec(MAXRES6)
3017       common /stochcalc/ stochforcvec
3018       logical lprn /.false./
3019       integer i,j
3020
3021 c      write (iout,*) "dc_old"
3022 c      do i=0,nres
3023 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
3024 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3025 c      enddo
3026       do j=1,3
3027         dc_work(j)=dc_old(j,0)
3028         d_t_work(j)=d_t_old(j,0)
3029         d_a_work(j)=d_a_old(j,0)
3030       enddo
3031       ind=3
3032       do i=nnt,nct-1
3033         do j=1,3
3034           dc_work(ind+j)=dc_old(j,i)
3035           d_t_work(ind+j)=d_t_old(j,i)
3036           d_a_work(ind+j)=d_a_old(j,i)
3037         enddo
3038         ind=ind+3
3039       enddo
3040       do i=nnt,nct
3041         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3042           do j=1,3
3043             dc_work(ind+j)=dc_old(j,i+nres)
3044             d_t_work(ind+j)=d_t_old(j,i+nres)
3045             d_a_work(ind+j)=d_a_old(j,i+nres)
3046           enddo
3047           ind=ind+3
3048         endif
3049       enddo
3050
3051 #ifndef LANG0
3052       if (lprn) then
3053       write (iout,*)
3054      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3055      &  " vrand_mat2"
3056       do i=1,dimen
3057         do j=1,dimen
3058           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3059      &      vfric_mat(i,j),afric_mat(i,j),
3060      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3061         enddo
3062       enddo
3063       endif
3064       do i=1,dimen
3065         ddt1=0.0d0
3066         ddt2=0.0d0
3067         do j=1,dimen
3068           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3069      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3070           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3071           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3072         enddo
3073         d_t_work_new(i)=ddt1+0.5d0*ddt2
3074         d_t_work(i)=ddt1+ddt2
3075       enddo
3076 #endif
3077       do j=1,3
3078         dc(j,0)=dc_work(j)
3079         d_t(j,0)=d_t_work(j)
3080       enddo
3081       ind=3
3082       do i=nnt,nct-1
3083         do j=1,3
3084           dc(j,i)=dc_work(ind+j)
3085           d_t(j,i)=d_t_work(ind+j)
3086         enddo
3087         ind=ind+3
3088       enddo
3089       do i=nnt,nct
3090         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3091           inres=i+nres
3092           do j=1,3
3093             dc(j,inres)=dc_work(ind+j)
3094             d_t(j,inres)=d_t_work(ind+j)
3095           enddo
3096           ind=ind+3
3097         endif
3098       enddo
3099       return
3100       end
3101 c--------------------------------------------------------------------------
3102       subroutine sd_verlet2_ciccotti
3103 c  Calculating the adjusted velocities for accelerations
3104       implicit none
3105       include 'DIMENSIONS'
3106       include 'COMMON.CONTROL'
3107       include 'COMMON.VAR'
3108       include 'COMMON.MD'
3109 #ifdef FIVEDIAG
3110       include 'COMMON.LAGRANGE.5diag'
3111 #else
3112       include 'COMMON.LAGRANGE'
3113 #endif
3114       include 'COMMON.LANGEVIN'
3115       include 'COMMON.CHAIN'
3116       include 'COMMON.DERIV'
3117       include 'COMMON.GEO'
3118       include 'COMMON.LOCAL'
3119       include 'COMMON.INTERACT'
3120       include 'COMMON.IOUNITS'
3121       include 'COMMON.NAMES'
3122       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3123       common /stochcalc/ stochforcvec
3124       integer i,j
3125 c
3126 c Compute the stochastic forces which contribute to velocity change
3127 c
3128       call stochastic_force(stochforcvecV)
3129 #ifndef LANG0
3130       do i=1,dimen
3131         ddt1=0.0d0
3132         ddt2=0.0d0
3133         do j=1,dimen
3134
3135           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3136 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3137           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3138         enddo
3139         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3140       enddo
3141 #endif
3142       do j=1,3
3143         d_t(j,0)=d_t_work(j)
3144       enddo
3145       ind=3
3146       do i=nnt,nct-1
3147         do j=1,3
3148           d_t(j,i)=d_t_work(ind+j)
3149         enddo
3150         ind=ind+3
3151       enddo
3152       do i=nnt,nct
3153         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3154           inres=i+nres
3155           do j=1,3
3156             d_t(j,inres)=d_t_work(ind+j)
3157           enddo
3158           ind=ind+3
3159         endif
3160       enddo
3161       return
3162       end
3163 #endif