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