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