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