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

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