source/unres/src_MD/old_F/MD.F

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