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               dc(:,0)=c(:,1)
1959               do i=1,nres-1
1960                 do j=1,3
1961                   dc(j,i)=c(j,i+1)-c(j,i)
1962                   dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1963                 enddo
1964               enddo
1965               do i=2,nres-1
1966                 do j=1,3
1967                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1968                   dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1969                 enddo
1970               enddo
1971               if (me.eq.king.or..not.out1file) then
1972               write (iout,*) "Energies before removing overlaps"
1973               call etotal(energia(0))
1974               call enerprint(energia(0))
1975               endif
1976 ! Remove SC overlaps if requested
1977               if (overlapsc) then
1978                 write (iout,*) 'Calling OVERLAP_SC'
1979                 call overlap_sc(fail)
1980                 if (fail) then
1981                   write (iout,*)
1982      &            "Failed to remove overlap from model",i_model
1983                   cycle
1984                 endif
1985               endif
1986               if (me.eq.king.or..not.out1file) then
1987               write (iout,*) "Energies after removing overlaps"
1988               call etotal(energia(0))
1989               call enerprint(energia(0))
1990               endif
1991 #ifdef SEARCHSC
1992 ! Search for better SC rotamers if requested
1993               if (searchsc) then
1994                 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1995                 print *,'SC_move',nft_sc,etot
1996                 if (me.eq.king.or..not.out1file)
1997      &            write(iout,*) 'SC_move',nft_sc,etot
1998               endif
1999               call etotal(energia(0))
2000 #endif
2001             enddo
2002             call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
2003      &        1,MPI_INTEGER,king,CG_COMM,IERROR)
2004             if (n_model_try.gt.nmodel_start .and.
2005      &         (me.eq.king .or. out1file)) then
2006               write (iout,*)
2007      &    "All models have irreparable overlaps. Trying randoms starts."
2008               iranconf=1
2009               i_model=nmodel_start+1
2010               goto 122
2011             endif
2012           else
2013 ! Remove SC overlaps if requested
2014               if (overlapsc) then
2015                 write (iout,*) 'Calling OVERLAP_SC'
2016                 call overlap_sc(fail)
2017                 if (fail) then
2018                   write (iout,*)
2019      &            "Failed to remove overlap"
2020                 endif
2021               endif
2022               if (me.eq.king.or..not.out1file) then
2023               write (iout,*) "Energies after removing overlaps"
2024               call etotal(energia(0))
2025               call enerprint(energia(0))
2026               endif
2027           endif
2028 C 8/22/17 AL Minimize initial structure
2029           if (dccart) then
2030             if (me.eq.king.or..not.out1file) write(iout,*)
2031      &        'Minimizing initial PDB structure: Calling MINIM_DC'
2032             call minim_dc(etot,iretcode,nfun)
2033           else
2034             call geom_to_var(nvar,varia)
2035             if(me.eq.king.or..not.out1file) write (iout,*)
2036      &        'Minimizing initial PDB structure: Calling MINIMIZE.'
2037             call minimize(etot,varia,iretcode,nfun)
2038             call var_to_geom(nvar,varia)
2039 #ifdef LBFGS
2040             if (me.eq.king.or..not.out1file)
2041      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2042             if(me.eq.king.or..not.out1file)
2043      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2044 #else
2045             if (me.eq.king.or..not.out1file)
2046      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2047             if(me.eq.king.or..not.out1file)
2048      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2049 #endif
2050           endif
2051         endif
2052         if (nmodel_start.gt.0 .and. me.eq.king) then
2053           write (iout,'(a)') "Task  Starting model"
2054           do i=0,nodes-1
2055             if (i_start_models(i).gt.nmodel_start) then
2056               write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
2057             else
2058               write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
2059      &          (:ilen(pdbfiles_chomo(i_start_models(i))))
2060             endif
2061           enddo
2062         endif
2063       endif ! .not. rest
2064       call chainbuild_cart
2065       call kinetic(EK)
2066       if (tbf) then
2067         call verlet_bath
2068       endif
2069       kinetic_T=2.0d0/(dimen3*Rb)*EK
2070       write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
2071       if(me.eq.king.or..not.out1file)then
2072        call cartprint
2073        call intout
2074       endif
2075 #ifdef MPI
2076       tt0=MPI_Wtime()
2077 #else
2078       tt0=tcpu()
2079 #endif
2080       call zerograd
2081       call etotal(potEcomp)
2082       if (large) call enerprint(potEcomp)
2083 #ifdef TIMING_ENE
2084 #ifdef MPI
2085       t_etotal=t_etotal+MPI_Wtime()-tt0
2086 #else
2087       t_etotal=t_etotal+tcpu()-tt0
2088 #endif
2089 #endif
2090       potE=potEcomp(0)
2091 c      write (iout,*) "PotE-homology",potE-potEcomp(27)
2092       call cartgrad
2093       call lagrangian
2094       call max_accel
2095       if (amax*d_time .gt. dvmax) then
2096         d_time=d_time*dvmax/amax
2097         if(me.eq.king.or..not.out1file) write (iout,*)
2098      &   "Time step reduced to",d_time,
2099      &   " because of too large initial acceleration."
2100       endif
2101       if(me.eq.king.or..not.out1file)then
2102        write(iout,*) "Potential energy and its components"
2103        call enerprint(potEcomp)
2104 c       write(iout,*) (potEcomp(i),i=0,n_ene)
2105       endif
2106       potE=potEcomp(0)-potEcomp(27)
2107 c      write (iout,*) "PotE-homology",potE
2108       totE=EK+potE
2109       itime=0
2110       itime_mat=itime
2111       if (ntwe.ne.0) call statout(itime)
2112       if(me.eq.king.or..not.out1file)
2113      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2114      &   " Kinetic energy",EK," Potential energy",potE,
2115      &   " Total energy",totE," Maximum acceleration ",
2116      &   amax
2117       if (large) then
2118         write (iout,*) "Initial coordinates"
2119         do i=1,nres
2120           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2121      &    (c(j,i+nres),j=1,3)
2122         enddo
2123         write (iout,*) "Initial dC"
2124         do i=0,nres
2125           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2126      &    (dc(j,i+nres),j=1,3)
2127         enddo
2128         write (iout,*) "Initial velocities"
2129         write (iout,"(13x,' backbone ',23x,' side chain')")
2130         do i=0,nres
2131           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2132      &    (d_t(j,i+nres),j=1,3)
2133         enddo
2134         write (iout,*) "Initial accelerations"
2135         do i=0,nres
2136 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2137           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2138      &    (d_a(j,i+nres),j=1,3)
2139         enddo
2140       endif
2141       do i=0,2*nres
2142         do j=1,3
2143           dc_old(j,i)=dc(j,i)
2144           d_t_old(j,i)=d_t(j,i)
2145           d_a_old(j,i)=d_a(j,i)
2146         enddo
2147 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2148       enddo
2149       if (RESPA) then
2150 #ifdef MPI
2151         tt0 =MPI_Wtime()
2152 #else
2153         tt0 = tcpu()
2154 #endif
2155         call zerograd
2156         call etotal_short(energia_short)
2157         if (large) call enerprint(potEcomp)
2158 #ifdef TIMING_ENE
2159 #ifdef MPI
2160         t_eshort=t_eshort+MPI_Wtime()-tt0
2161 #else
2162         t_eshort=t_eshort+tcpu()-tt0
2163 #endif
2164 #endif
2165         call cartgrad
2166         call lagrangian
2167         if(.not.out1file .and. large) then
2168           write (iout,*) "energia_long",energia_long(0),
2169      &     " energia_short",energia_short(0),
2170      &     " total",energia_long(0)+energia_short(0)
2171           write (iout,*) "Initial fast-force accelerations"
2172           do i=0,nres
2173             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2174      &      (d_a(j,i+nres),j=1,3)
2175           enddo
2176         endif
2177 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2178         do i=0,2*nres
2179           do j=1,3
2180             d_a_short(j,i)=d_a(j,i)
2181           enddo
2182         enddo
2183 #ifdef MPI
2184         tt0=MPI_Wtime()
2185 #else
2186         tt0=tcpu()
2187 #endif
2188         call zerograd
2189         call etotal_long(energia_long)
2190         if (large) call enerprint(potEcomp)
2191 #ifdef TIMING_ENE
2192 #ifdef MPI
2193         t_elong=t_elong+MPI_Wtime()-tt0
2194 #else
2195         t_elong=t_elong+tcpu()-tt0
2196 #endif
2197 #endif
2198         call cartgrad
2199         call lagrangian
2200         if(.not.out1file .and. large) then
2201           write (iout,*) "energia_long",energia_long(0)
2202           write (iout,*) "Initial slow-force accelerations"
2203           do i=0,nres
2204             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2205      &      (d_a(j,i+nres),j=1,3)
2206           enddo
2207         endif
2208 #ifdef MPI
2209         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2210 #else
2211         t_enegrad=t_enegrad+tcpu()-tt0
2212 #endif
2213       endif
2214       return
2215       end
2216 c-----------------------------------------------------------
2217       subroutine random_vel
2218       implicit none
2219       include 'DIMENSIONS'
2220       include 'COMMON.CONTROL'
2221       include 'COMMON.VAR'
2222       include 'COMMON.MD'
2223 #ifdef FIVEDIAG
2224        include 'COMMON.LAGRANGE.5diag'
2225 #else
2226        include 'COMMON.LAGRANGE'
2227 #endif
2228 #ifndef LANG0
2229       include 'COMMON.LANGEVIN'
2230 #else
2231 #ifdef FIVEDIAG
2232       include 'COMMON.LANGEVIN.lang0.5diag'
2233 #else
2234       include 'COMMON.LANGEVIN.lang0'
2235 #endif
2236 #endif
2237       include 'COMMON.CHAIN'
2238       include 'COMMON.DERIV'
2239       include 'COMMON.GEO'
2240       include 'COMMON.LOCAL'
2241       include 'COMMON.INTERACT'
2242       include 'COMMON.IOUNITS'
2243       include 'COMMON.NAMES'
2244       include 'COMMON.TIME1'
2245       double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2246       integer i,ii,j,k,l,ind
2247       double precision anorm_distr
2248       logical lprn /.false./
2249 #ifdef FIVEDIAG
2250       integer ichain,n,innt,inct,ibeg,ierr
2251       double precision work(8*maxres6)
2252       integer iwork(maxres6)
2253       double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2254      & Gvec(maxres2_chain,maxres2_chain)
2255       common /przechowalnia/Ghalf,Geigen,Gvec
2256 #ifdef DEBUG
2257       double precision inertia(maxres2_chain,maxres2_chain)
2258 #endif
2259 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2260 c First generate velocities in the eigenspace of the G matrix
2261 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2262 c      call flush(iout)
2263 #ifdef DEBUG
2264       write (iout,*) "Random_vel, fivediag"
2265 #endif
2266       d_t=0.0d0
2267       Ek2=0.0d0
2268       EK=0.0d0
2269       Ek3=0.0d0
2270       do ichain=1,nchain
2271         ind=0
2272         ghalf=0.0d0
2273         n=dimen_chain(ichain)
2274         innt=iposd_chain(ichain)
2275         inct=innt+n-1
2276 #ifdef DEBUG
2277         write (iout,*) "Chain",ichain," n",n," start",innt
2278         do i=innt,inct
2279           if (i.lt.inct-1) then
2280             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2281      &         DU2orig(i)
2282           else if (i.eq.inct-1) then
2283             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2284           else
2285             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2286           endif
2287         enddo
2288 #endif
2289         ghalf(ind+1)=dmorig(innt)
2290         ghalf(ind+2)=du1orig(innt)
2291         ghalf(ind+3)=dmorig(innt+1)
2292         ind=ind+3
2293         do i=3,n
2294           ind=ind+i-3
2295 c          write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2296 c     &       " indu1",innt+i-1," indm",innt+i
2297           ghalf(ind+1)=du2orig(innt-1+i-2)
2298           ghalf(ind+2)=du1orig(innt-1+i-1)
2299           ghalf(ind+3)=dmorig(innt-1+i)
2300 c          write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2301 c     &       "DU1",innt-1+i-1,"DM ",innt-1+i
2302           ind=ind+3
2303         enddo
2304 #ifdef DEBUG
2305         ind=0
2306         do i=1,n
2307           do j=1,i
2308             ind=ind+1
2309             inertia(i,j)=ghalf(ind)
2310             inertia(j,i)=ghalf(ind)
2311           enddo
2312         enddo
2313 #endif
2314 #ifdef DEBUG
2315         write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2316         write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2317         call matoutr(n,ghalf)
2318 #endif
2319         call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2320         if (large) then
2321           write (iout,'(//a,i3)')
2322      &    "Eigenvectors and eigenvalues of the G matrix chain",ichain
2323           call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2324         endif
2325 #ifdef DIAGCHECK
2326 c check diagonalization
2327         do i=1,n
2328           do j=1,n
2329             aux=0.0d0
2330             do k=1,n
2331               do l=1,n
2332                 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2333               enddo
2334             enddo
2335             if (i.eq.j) then
2336               write (iout,*) i,j,aux,geigen(i)
2337             else
2338               write (iout,*) i,j,aux
2339             endif
2340           enddo
2341         enddo
2342 #endif
2343         xv=0.0d0
2344         ii=0
2345         do i=1,n
2346           do k=1,3
2347             ii=ii+1
2348             sigv=dsqrt((Rb*t_bath)/geigen(i))
2349             lowb=-5*sigv
2350             highb=5*sigv
2351             d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2352             EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2353 c            write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2354 c     &      " d_t_work_new",d_t_work_new(ii)
2355           enddo
2356         enddo
2357         do k=1,3
2358           do i=1,n
2359             ind=(i-1)*3+k
2360             d_t_work(ind)=0.0d0
2361             do j=1,n
2362               d_t_work(ind)=d_t_work(ind)
2363      &                +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2364             enddo
2365 c            write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2366 c            call flush(iout)
2367           enddo
2368         enddo
2369 #ifdef DEBUG
2370         aux=0.0d0
2371         do k=1,3
2372           do i=1,n
2373             do j=1,n
2374             aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2375             enddo
2376           enddo
2377         enddo
2378         Ek3=Ek3+aux/2
2379 #endif
2380 c Transfer to the d_t vector
2381         innt=chain_border(1,ichain)
2382         inct=chain_border(2,ichain)
2383         ind=0
2384 c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
2385         do i=innt,inct
2386           do j=1,3
2387             ind=ind+1
2388             d_t(j,i)=d_t_work(ind)
2389           enddo
2390           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2391             do j=1,3
2392               ind=ind+1
2393               d_t(j,i+nres)=d_t_work(ind)
2394             enddo
2395           endif
2396         enddo
2397       enddo
2398       if (large) then
2399         write (iout,*)
2400         write (iout,*) "Random velocities in the Calpha,SC space"
2401         do i=1,nres
2402           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2403      &    restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2404         enddo
2405       endif
2406       call kinetic_CASC(Ek1)
2407 !
2408 ! Transform the velocities to virtual-bond space
2409 !
2410 #define WLOS
2411 #ifdef WLOS
2412       if (nnt.eq.1) then
2413         d_t(:,0)=d_t(:,1)
2414       endif
2415       do i=1,nres
2416         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2417           do j=1,3
2418             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2419           enddo
2420         else
2421           do j=1,3
2422             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2423             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2424           enddo
2425         end if
2426       enddo
2427       d_t(:,nres)=0.0d0
2428       d_t(:,nct)=0.0d0
2429       d_t(:,2*nres)=0.0d0
2430       if (nnt.gt.1) then
2431         d_t(:,0)=d_t(:,1)
2432         d_t(:,1)=0.0d0
2433       endif
2434 c      d_a(:,0)=d_a(:,1)
2435 c      d_a(:,1)=0.0d0
2436 c      write (iout,*) "Shifting accelerations"
2437       do ichain=2,nchain
2438 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
2439 c     &     chain_border1(1,ichain)
2440         d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2441         d_t(:,chain_border1(1,ichain))=0.0d0
2442       enddo
2443 c      write (iout,*) "Adding accelerations"
2444       do ichain=2,nchain
2445 c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2446 c     &   chain_border(2,ichain-1)
2447         d_t(:,chain_border1(1,ichain)-1)=
2448      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2449         d_t(:,chain_border(2,ichain-1))=0.0d0
2450       enddo
2451       do ichain=2,nchain
2452         write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2453      &   chain_border(2,ichain-1)
2454         d_t(:,chain_border1(1,ichain)-1)=
2455      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2456         d_t(:,chain_border(2,ichain-1))=0.0d0
2457       enddo
2458 #else
2459       ibeg=0
2460 c      do j=1,3
2461 c        d_t(j,0)=d_t(j,nnt)
2462 c      enddo
2463       do ichain=1,nchain
2464       innt=chain_border(1,ichain)
2465       inct=chain_border(2,ichain)
2466 c      write (iout,*) "ichain",ichain," innt",innt," inct",inct
2467 c      write (iout,*) "ibeg",ibeg
2468       do j=1,3
2469         d_t(j,ibeg)=d_t(j,innt)
2470       enddo
2471       ibeg=inct+1
2472       do i=innt,inct
2473         if (iabs(itype(i).eq.10)) then
2474 c          write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2475           do j=1,3
2476             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2477           enddo
2478         else
2479           do j=1,3
2480             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2481             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2482           enddo
2483         end if
2484       enddo
2485       enddo
2486 #endif
2487       if (large) then
2488         write (iout,*)
2489         write (iout,*)
2490      &    "Random velocities in the virtual-bond-vector space"
2491         write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2492         do i=1,nres
2493           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2494      &     restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2495         enddo
2496         write (iout,*)
2497         write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2498      &   Ek
2499         write (iout,*)
2500      &   "Kinetic temperatures from inertia matrix eigenvalues",
2501      &   2*Ek/(3*dimen*Rb)
2502 #ifdef DEBUG
2503         write (iout,*) "Kinetic energy from inertia matrix",Ek3
2504         write (iout,*) "Kinetic temperatures from inertia",
2505      &   2*Ek3/(3*dimen*Rb)
2506 #endif
2507         write (iout,*) "Kinetic energy from velocities in CA-SC space",
2508      &   Ek1
2509         write (iout,*)
2510      &   "Kinetic temperatures from velovities in CA-SC space",
2511      &   2*Ek1/(3*dimen*Rb)
2512         call kinetic(Ek1)
2513         write (iout,*)
2514      &   "Kinetic energy from virtual-bond-vector velocities",Ek1
2515         write (iout,*)
2516      &   "Kinetic temperature from virtual-bond-vector velocities ",
2517      &   2*Ek1/(dimen3*Rb)
2518       endif
2519 #else
2520       xv=0.0d0
2521       ii=0
2522       do i=1,dimen
2523         do k=1,3
2524           ii=ii+1
2525           sigv=dsqrt((Rb*t_bath)/geigen(i))
2526           lowb=-5*sigv
2527           highb=5*sigv
2528           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2529
2530 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2531 c     &      " d_t_work_new",d_t_work_new(ii)
2532         enddo
2533       enddo
2534 C       if (SELFGUIDE.gt.0) then
2535 C       distance=0.0
2536 C       do j=1,3
2537 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2538 C       distance=distance+vec_afm(j)**2
2539 C       enddo
2540 C       distance=dsqrt(distance)
2541 C       do j=1,3
2542 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2543 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2544 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2545 C     &    d_t_work_new(j+(afmend-1)*3)
2546 C       enddo
2547
2548 C       endif
2549
2550 c diagnostics
2551 c      Ek1=0.0d0
2552 c      ii=0
2553 c      do i=1,dimen
2554 c        do k=1,3
2555 c          ii=ii+1
2556 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2557 c        enddo
2558 c      enddo
2559 c      write (iout,*) "Ek from eigenvectors",Ek1
2560 c end diagnostics
2561 c Transform velocities to UNRES coordinate space
2562       do k=0,2
2563         do i=1,dimen
2564           ind=(i-1)*3+k+1
2565           d_t_work(ind)=0.0d0
2566           do j=1,dimen
2567             d_t_work(ind)=d_t_work(ind)
2568      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2569           enddo
2570 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2571 c          call flush(iout)
2572         enddo
2573       enddo
2574 c Transfer to the d_t vector
2575       do j=1,3
2576         d_t(j,0)=d_t_work(j)
2577       enddo
2578       ind=3
2579       do i=nnt,nct-1
2580         do j=1,3
2581           ind=ind+1
2582           d_t(j,i)=d_t_work(ind)
2583         enddo
2584       enddo
2585       do i=nnt,nct
2586         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2587           do j=1,3
2588             ind=ind+1
2589             d_t(j,i+nres)=d_t_work(ind)
2590           enddo
2591         endif
2592       enddo
2593 c      call kinetic(EK)
2594 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2595 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2596 c      call flush(iout)
2597 #endif
2598       return
2599       end
2600 #ifndef LANG0
2601 c-----------------------------------------------------------
2602       subroutine sd_verlet_p_setup
2603 c Sets up the parameters of stochastic Verlet algorithm
2604       implicit none
2605       include 'DIMENSIONS'
2606 #ifdef MPI
2607       include 'mpif.h'
2608 #endif
2609       include 'COMMON.CONTROL'
2610       include 'COMMON.VAR'
2611       include 'COMMON.MD'
2612 #ifndef LANG0
2613       include 'COMMON.LANGEVIN'
2614 #else
2615 #ifdef FIVEDIAG
2616       include 'COMMON.LANGEVIN.lang0.5diag'
2617 #else
2618       include 'COMMON.LANGEVIN.lang0'
2619 #endif
2620 #endif
2621       include 'COMMON.CHAIN'
2622       include 'COMMON.DERIV'
2623       include 'COMMON.GEO'
2624       include 'COMMON.LOCAL'
2625       include 'COMMON.INTERACT'
2626       include 'COMMON.IOUNITS'
2627       include 'COMMON.NAMES'
2628       include 'COMMON.TIME1'
2629       double precision emgdt(MAXRES6),
2630      & pterm,vterm,rho,rhoc,vsig,
2631      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2632      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2633      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2634       logical lprn /.false./
2635       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2636       double precision ktm
2637 #ifdef MPI
2638       tt0 = MPI_Wtime()
2639 #else
2640       tt0 = tcpu()
2641 #endif
2642 c
2643 c AL 8/17/04 Code adapted from tinker
2644 c
2645 c Get the frictional and random terms for stochastic dynamics in the
2646 c eigenspace of mass-scaled UNRES friction matrix
2647 c
2648       do i = 1, dimen
2649             gdt = fricgam(i) * d_time
2650 c
2651 c Stochastic dynamics reduces to simple MD for zero friction
2652 c
2653             if (gdt .le. zero) then
2654                pfric_vec(i) = 1.0d0
2655                vfric_vec(i) = d_time
2656                afric_vec(i) = 0.5d0 * d_time * d_time
2657                prand_vec(i) = 0.0d0
2658                vrand_vec1(i) = 0.0d0
2659                vrand_vec2(i) = 0.0d0
2660 c
2661 c Analytical expressions when friction coefficient is large
2662 c
2663             else
2664                if (gdt .ge. gdt_radius) then
2665                   egdt = dexp(-gdt)
2666                   pfric_vec(i) = egdt
2667                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2668                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2669                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2670                   vterm = 1.0d0 - egdt**2
2671                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2672 c
2673 c Use series expansions when friction coefficient is small
2674 c
2675                else
2676                   gdt2 = gdt * gdt
2677                   gdt3 = gdt * gdt2
2678                   gdt4 = gdt2 * gdt2
2679                   gdt5 = gdt2 * gdt3
2680                   gdt6 = gdt3 * gdt3
2681                   gdt7 = gdt3 * gdt4
2682                   gdt8 = gdt4 * gdt4
2683                   gdt9 = gdt4 * gdt5
2684                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2685      &                          - gdt5/120.0d0 + gdt6/720.0d0
2686      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2687      &                          - gdt9/362880.0d0) / fricgam(i)**2
2688                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2689                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2690                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2691      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2692      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2693      &                       + 127.0d0*gdt9/90720.0d0
2694                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2695      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2696      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2697      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2698                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2699      &                       - 17.0d0*gdt2/1280.0d0
2700      &                       + 17.0d0*gdt3/6144.0d0
2701      &                       + 40967.0d0*gdt4/34406400.0d0
2702      &                       - 57203.0d0*gdt5/275251200.0d0
2703      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2704                end if
2705 c
2706 c Compute the scaling factors of random terms for the nonzero friction case
2707 c
2708                ktm = 0.5d0*d_time/fricgam(i)
2709                psig = dsqrt(ktm*pterm) / fricgam(i)
2710                vsig = dsqrt(ktm*vterm)
2711                rhoc = dsqrt(1.0d0 - rho*rho)
2712                prand_vec(i) = psig
2713                vrand_vec1(i) = vsig * rho
2714                vrand_vec2(i) = vsig * rhoc
2715             end if
2716       end do
2717       if (lprn) then
2718       write (iout,*)
2719      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2720      &  " vrand_vec2"
2721       do i=1,dimen
2722         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2723      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2724       enddo
2725       endif
2726 c
2727 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2728 c
2729       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2730       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2731       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2732       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2733       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2734       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2735 #ifdef MPI
2736       t_sdsetup=t_sdsetup+MPI_Wtime()
2737 #else
2738       t_sdsetup=t_sdsetup+tcpu()-tt0
2739 #endif
2740       return
2741       end
2742 c-------------------------------------------------------------
2743       subroutine eigtransf1(n,ndim,ab,d,c)
2744       implicit none
2745       integer n,ndim
2746       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2747       integer i,j,k
2748       do i=1,n
2749         do j=1,n
2750           c(i,j)=0.0d0
2751           do k=1,n
2752             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2753           enddo
2754         enddo
2755       enddo
2756       return
2757       end
2758 c-------------------------------------------------------------
2759       subroutine eigtransf(n,ndim,a,b,d,c)
2760       implicit none
2761       integer n,ndim
2762       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2763       integer i,j,k
2764       do i=1,n
2765         do j=1,n
2766           c(i,j)=0.0d0
2767           do k=1,n
2768             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2769           enddo
2770         enddo
2771       enddo
2772       return
2773       end
2774 c-------------------------------------------------------------
2775       subroutine sd_verlet1
2776 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2777       implicit none
2778       include 'DIMENSIONS'
2779       include 'COMMON.CONTROL'
2780       include 'COMMON.VAR'
2781       include 'COMMON.MD'
2782 #ifdef FIVEDIAG
2783        include 'COMMON.LAGRANGE.5diag'
2784 #else
2785        include 'COMMON.LAGRANGE'
2786 #endif
2787 #ifndef LANG0
2788       include 'COMMON.LANGEVIN'
2789 #else
2790 #ifdef FIVEDIAG
2791       include 'COMMON.LANGEVIN.lang0.5diag'
2792 #else
2793       include 'COMMON.LANGEVIN.lang0'
2794 #endif
2795 #endif
2796       include 'COMMON.CHAIN'
2797       include 'COMMON.DERIV'
2798       include 'COMMON.GEO'
2799       include 'COMMON.LOCAL'
2800       include 'COMMON.INTERACT'
2801       include 'COMMON.IOUNITS'
2802       include 'COMMON.NAMES'
2803       double precision stochforcvec(MAXRES6)
2804       common /stochcalc/ stochforcvec
2805       logical lprn /.false./
2806       integer i,j,ind,inres
2807
2808 c      write (iout,*) "dc_old"
2809 c      do i=0,nres
2810 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2811 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2812 c      enddo
2813       do j=1,3
2814         dc_work(j)=dc_old(j,0)
2815         d_t_work(j)=d_t_old(j,0)
2816         d_a_work(j)=d_a_old(j,0)
2817       enddo
2818       ind=3
2819       do i=nnt,nct-1
2820         do j=1,3
2821           dc_work(ind+j)=dc_old(j,i)
2822           d_t_work(ind+j)=d_t_old(j,i)
2823           d_a_work(ind+j)=d_a_old(j,i)
2824         enddo
2825         ind=ind+3
2826       enddo
2827       do i=nnt,nct
2828         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2829           do j=1,3
2830             dc_work(ind+j)=dc_old(j,i+nres)
2831             d_t_work(ind+j)=d_t_old(j,i+nres)
2832             d_a_work(ind+j)=d_a_old(j,i+nres)
2833           enddo
2834           ind=ind+3
2835         endif
2836       enddo
2837 #ifndef LANG0
2838       if (lprn) then
2839       write (iout,*)
2840      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2841      &  " vrand_mat2"
2842       do i=1,dimen
2843         do j=1,dimen
2844           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2845      &      vfric_mat(i,j),afric_mat(i,j),
2846      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2847         enddo
2848       enddo
2849       endif
2850       do i=1,dimen
2851         ddt1=0.0d0
2852         ddt2=0.0d0
2853         do j=1,dimen
2854           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2855      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2856           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2857           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2858         enddo
2859         d_t_work_new(i)=ddt1+0.5d0*ddt2
2860         d_t_work(i)=ddt1+ddt2
2861       enddo
2862 #endif
2863       do j=1,3
2864         dc(j,0)=dc_work(j)
2865         d_t(j,0)=d_t_work(j)
2866       enddo
2867       ind=3
2868       do i=nnt,nct-1
2869         do j=1,3
2870           dc(j,i)=dc_work(ind+j)
2871           d_t(j,i)=d_t_work(ind+j)
2872         enddo
2873         ind=ind+3
2874       enddo
2875       do i=nnt,nct
2876         if (itype(i).ne.10) then
2877           inres=i+nres
2878           do j=1,3
2879             dc(j,inres)=dc_work(ind+j)
2880             d_t(j,inres)=d_t_work(ind+j)
2881           enddo
2882           ind=ind+3
2883         endif
2884       enddo
2885       return
2886       end
2887 c--------------------------------------------------------------------------
2888       subroutine sd_verlet2
2889 c  Calculating the adjusted velocities for accelerations
2890       implicit none
2891       include 'DIMENSIONS'
2892       include 'COMMON.CONTROL'
2893       include 'COMMON.VAR'
2894       include 'COMMON.MD'
2895 #ifdef FIVEDIAG
2896        include 'COMMON.LAGRANGE.5diag'
2897 #else
2898        include 'COMMON.LAGRANGE'
2899 #endif
2900 #ifndef LANG0
2901       include 'COMMON.LANGEVIN'
2902 #else
2903 #ifdef FIVEDIAG
2904       include 'COMMON.LANGEVIN.lang0.5diag'
2905 #else
2906       include 'COMMON.LANGEVIN.lang0'
2907 #endif
2908 #endif
2909       include 'COMMON.CHAIN'
2910       include 'COMMON.DERIV'
2911       include 'COMMON.GEO'
2912       include 'COMMON.LOCAL'
2913       include 'COMMON.INTERACT'
2914       include 'COMMON.IOUNITS'
2915       include 'COMMON.NAMES'
2916       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2917       common /stochcalc/ stochforcvec
2918       integer i,j,ind,inres
2919 c
2920 c Compute the stochastic forces which contribute to velocity change
2921 c
2922       call stochastic_force(stochforcvecV)
2923
2924 #ifndef LANG0
2925       do i=1,dimen
2926         ddt1=0.0d0
2927         ddt2=0.0d0
2928         do j=1,dimen
2929           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2930           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2931      &     vrand_mat2(i,j)*stochforcvecV(j)
2932         enddo
2933         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2934       enddo
2935 #endif
2936       do j=1,3
2937         d_t(j,0)=d_t_work(j)
2938       enddo
2939       ind=3
2940       do i=nnt,nct-1
2941         do j=1,3
2942           d_t(j,i)=d_t_work(ind+j)
2943         enddo
2944         ind=ind+3
2945       enddo
2946       do i=nnt,nct
2947         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2948           inres=i+nres
2949           do j=1,3
2950             d_t(j,inres)=d_t_work(ind+j)
2951           enddo
2952           ind=ind+3
2953         endif
2954       enddo
2955       return
2956       end
2957 c-----------------------------------------------------------
2958       subroutine sd_verlet_ciccotti_setup
2959 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2960 c version
2961       implicit none
2962       include 'DIMENSIONS'
2963 #ifdef MPI
2964       include 'mpif.h'
2965 #endif
2966       include 'COMMON.CONTROL'
2967       include 'COMMON.VAR'
2968       include 'COMMON.MD'
2969 #ifdef FIVEDIAG
2970       include 'COMMON.LAGRANGE.5diag'
2971 #else
2972       include 'COMMON.LAGRANGE'
2973 #endif
2974       include 'COMMON.LANGEVIN'
2975       include 'COMMON.CHAIN'
2976       include 'COMMON.DERIV'
2977       include 'COMMON.GEO'
2978       include 'COMMON.LOCAL'
2979       include 'COMMON.INTERACT'
2980       include 'COMMON.IOUNITS'
2981       include 'COMMON.NAMES'
2982       include 'COMMON.TIME1'
2983       double precision emgdt(MAXRES6),
2984      & pterm,vterm,rho,rhoc,vsig,
2985      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2986      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2987      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2988       logical lprn /.false./
2989       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2990       double precision ktm
2991       integer i
2992 #ifdef MPI
2993       tt0 = MPI_Wtime()
2994 #else
2995       tt0 = tcpu()
2996 #endif
2997 c
2998 c AL 8/17/04 Code adapted from tinker
2999 c
3000 c Get the frictional and random terms for stochastic dynamics in the
3001 c eigenspace of mass-scaled UNRES friction matrix
3002 c
3003       do i = 1, dimen
3004             write (iout,*) "i",i," fricgam",fricgam(i)
3005             gdt = fricgam(i) * d_time
3006 c
3007 c Stochastic dynamics reduces to simple MD for zero friction
3008 c
3009             if (gdt .le. zero) then
3010                pfric_vec(i) = 1.0d0
3011                vfric_vec(i) = d_time
3012                afric_vec(i) = 0.5d0*d_time*d_time
3013                prand_vec(i) = afric_vec(i)
3014                vrand_vec2(i) = vfric_vec(i)
3015 c
3016 c Analytical expressions when friction coefficient is large
3017 c
3018             else
3019                egdt = dexp(-gdt)
3020                pfric_vec(i) = egdt
3021                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3022                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3023                prand_vec(i) = afric_vec(i)
3024                vrand_vec2(i) = vfric_vec(i)
3025 c
3026 c Compute the scaling factors of random terms for the nonzero friction case
3027 c
3028 c               ktm = 0.5d0*d_time/fricgam(i)
3029 c               psig = dsqrt(ktm*pterm) / fricgam(i)
3030 c               vsig = dsqrt(ktm*vterm)
3031 c               prand_vec(i) = psig*afric_vec(i)
3032 c               vrand_vec2(i) = vsig*vfric_vec(i)
3033             end if
3034       end do
3035       if (lprn) then
3036       write (iout,*)
3037      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
3038      &  " vrand_vec2"
3039       do i=1,dimen
3040         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
3041      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3042       enddo
3043       endif
3044 c
3045 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3046 c
3047       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3048       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3049       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3050       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3051       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3052 #ifdef MPI
3053       t_sdsetup=t_sdsetup+MPI_Wtime()
3054 #else
3055       t_sdsetup=t_sdsetup+tcpu()-tt0
3056 #endif
3057       return
3058       end
3059 c-------------------------------------------------------------
3060       subroutine sd_verlet1_ciccotti
3061 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3062       implicit none
3063       include 'DIMENSIONS'
3064 #ifdef MPI
3065       include 'mpif.h'
3066 #endif
3067       include 'COMMON.CONTROL'
3068       include 'COMMON.VAR'
3069       include 'COMMON.MD'
3070 #ifdef FIVEDIAG
3071       include 'COMMON.LAGRANGE.5diag'
3072 #else
3073       include 'COMMON.LAGRANGE'
3074 #endif
3075       include 'COMMON.LANGEVIN'
3076       include 'COMMON.CHAIN'
3077       include 'COMMON.DERIV'
3078       include 'COMMON.GEO'
3079       include 'COMMON.LOCAL'
3080       include 'COMMON.INTERACT'
3081       include 'COMMON.IOUNITS'
3082       include 'COMMON.NAMES'
3083       double precision stochforcvec(MAXRES6)
3084       common /stochcalc/ stochforcvec
3085       logical lprn /.false./
3086       integer i,j
3087
3088 c      write (iout,*) "dc_old"
3089 c      do i=0,nres
3090 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
3091 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3092 c      enddo
3093       do j=1,3
3094         dc_work(j)=dc_old(j,0)
3095         d_t_work(j)=d_t_old(j,0)
3096         d_a_work(j)=d_a_old(j,0)
3097       enddo
3098       ind=3
3099       do i=nnt,nct-1
3100         do j=1,3
3101           dc_work(ind+j)=dc_old(j,i)
3102           d_t_work(ind+j)=d_t_old(j,i)
3103           d_a_work(ind+j)=d_a_old(j,i)
3104         enddo
3105         ind=ind+3
3106       enddo
3107       do i=nnt,nct
3108         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3109           do j=1,3
3110             dc_work(ind+j)=dc_old(j,i+nres)
3111             d_t_work(ind+j)=d_t_old(j,i+nres)
3112             d_a_work(ind+j)=d_a_old(j,i+nres)
3113           enddo
3114           ind=ind+3
3115         endif
3116       enddo
3117
3118 #ifndef LANG0
3119       if (lprn) then
3120       write (iout,*)
3121      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3122      &  " vrand_mat2"
3123       do i=1,dimen
3124         do j=1,dimen
3125           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3126      &      vfric_mat(i,j),afric_mat(i,j),
3127      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3128         enddo
3129       enddo
3130       endif
3131       do i=1,dimen
3132         ddt1=0.0d0
3133         ddt2=0.0d0
3134         do j=1,dimen
3135           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3136      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3137           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3138           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3139         enddo
3140         d_t_work_new(i)=ddt1+0.5d0*ddt2
3141         d_t_work(i)=ddt1+ddt2
3142       enddo
3143 #endif
3144       do j=1,3
3145         dc(j,0)=dc_work(j)
3146         d_t(j,0)=d_t_work(j)
3147       enddo
3148       ind=3
3149       do i=nnt,nct-1
3150         do j=1,3
3151           dc(j,i)=dc_work(ind+j)
3152           d_t(j,i)=d_t_work(ind+j)
3153         enddo
3154         ind=ind+3
3155       enddo
3156       do i=nnt,nct
3157         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3158           inres=i+nres
3159           do j=1,3
3160             dc(j,inres)=dc_work(ind+j)
3161             d_t(j,inres)=d_t_work(ind+j)
3162           enddo
3163           ind=ind+3
3164         endif
3165       enddo
3166       return
3167       end
3168 c--------------------------------------------------------------------------
3169       subroutine sd_verlet2_ciccotti
3170 c  Calculating the adjusted velocities for accelerations
3171       implicit none
3172       include 'DIMENSIONS'
3173       include 'COMMON.CONTROL'
3174       include 'COMMON.VAR'
3175       include 'COMMON.MD'
3176 #ifdef FIVEDIAG
3177       include 'COMMON.LAGRANGE.5diag'
3178 #else
3179       include 'COMMON.LAGRANGE'
3180 #endif
3181       include 'COMMON.LANGEVIN'
3182       include 'COMMON.CHAIN'
3183       include 'COMMON.DERIV'
3184       include 'COMMON.GEO'
3185       include 'COMMON.LOCAL'
3186       include 'COMMON.INTERACT'
3187       include 'COMMON.IOUNITS'
3188       include 'COMMON.NAMES'
3189       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3190       common /stochcalc/ stochforcvec
3191       integer i,j
3192 c
3193 c Compute the stochastic forces which contribute to velocity change
3194 c
3195       call stochastic_force(stochforcvecV)
3196 #ifndef LANG0
3197       do i=1,dimen
3198         ddt1=0.0d0
3199         ddt2=0.0d0
3200         do j=1,dimen
3201
3202           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3203 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3204           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3205         enddo
3206         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3207       enddo
3208 #endif
3209       do j=1,3
3210         d_t(j,0)=d_t_work(j)
3211       enddo
3212       ind=3
3213       do i=nnt,nct-1
3214         do j=1,3
3215           d_t(j,i)=d_t_work(ind+j)
3216         enddo
3217         ind=ind+3
3218       enddo
3219       do i=nnt,nct
3220         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3221           inres=i+nres
3222           do j=1,3
3223             d_t(j,inres)=d_t_work(ind+j)
3224           enddo
3225           ind=ind+3
3226         endif
3227       enddo
3228       return
3229       end
3230 #endif