source/unres/MD.f90

   1       module MDyn
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use math
   6       use md_calc
   7       use geometry_data
   8       use io_base
   9       use geometry
  10       use energy
  11       use MD_data
  12       use REMD
  13
  14       implicit none
  15 !-----------------------------------------------------------------------------
  16 ! common.MD
  17 !      common /mdgrad/ in module.energy
  18 !      common /back_constr/ in module.energy
  19 !      common /qmeas/ in module.energy
  20 !      common /mdpar/
  21 !      common /MDcalc/
  22 !      common /lagrange/
  23       real(kind=8),dimension(:),allocatable :: d_t_work,&
  24        d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
  25       real(kind=8),dimension(:,:),allocatable :: d_t_new,&
  26        d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
  27 !      real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
  28 !      real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
  29 !       Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
  30 !      real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
  31 !      integer :: dimen,dimen1,dimen3
  32 !      integer :: lang,count_reset_moment,count_reset_vel
  33 !      logical :: reset_moment,reset_vel,rattle,RESPA
  34 !      common /inertia/
  35 !      common /langevin/
  36 !      real(kind=8) :: rwat,etawat,stdfp,cPoise
  37 !      real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
  38 !      real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
  39       real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
  40 !-----------------------------------------------------------------------------
  41 ! 'sizes.i'
  42 !
  43 !
  44 !     ###################################################
  45 !     ##  COPYRIGHT (C)  1992  by  Jay William Ponder  ##
  46 !     ##              All Rights Reserved              ##
  47 !     ###################################################
  48 !
  49 !     #############################################################
  50 !     ##                                                         ##
  51 !     ##  sizes.i  --  parameter values to set array dimensions  ##
  52 !     ##                                                         ##
  53 !     #############################################################
  54 !
  55 !
  56 !     "sizes.i" sets values for critical array dimensions used
  57 !     throughout the software; these parameters will fix the size
  58 !     of the largest systems that can be handled; values too large
  59 !     for the computer's memory and/or swap space to accomodate
  60 !     will result in poor performance or outright failure
  61 !
  62 !     parameter:      maximum allowed number of:
  63 !
  64 !     maxatm          atoms in the molecular system
  65 !     maxval          atoms directly bonded to an atom
  66 !     maxgrp       !  user-defined groups of atoms
  67 !     maxtyp          force field atom type definitions
  68 !     maxclass        force field atom class definitions
  69 !     maxkey          lines in the keyword file
  70 !     maxrot          bonds for torsional rotation
  71 !     maxvar          optimization variables (vector storage)
  72 !     maxopt          optimization variables (matrix storage)
  73 !     maxhess         off-diagonal Hessian elements
  74 !     maxlight        sites for method of lights neighbors
  75 !     maxvib          vibrational frequencies
  76 !     maxgeo          distance geometry points
  77 !     maxcell         unit cells in replicated crystal
  78 !     maxring         3-, 4-, or 5-membered rings
  79 !     maxfix          geometric restraints
  80 !     maxbio          biopolymer atom definitions
  81 !     maxres          residues in the macromolecule
  82 !     maxamino        amino acid residue types
  83 !     maxnuc          nucleic acid residue types
  84 !     maxbnd          covalent bonds in molecular system
  85 !     maxang          bond angles in molecular system
  86 !     maxtors         torsional angles in molecular system
  87 !     maxpi           atoms in conjugated pisystem
  88 !     maxpib          covalent bonds involving pisystem
  89 !     maxpit          torsional angles involving pisystem
  90 !
  91 !
  92 !el      integer maxatm,maxval,maxgrp
  93 !el      integer maxtyp,maxclass,maxkey
  94 !el      integer maxrot,maxopt
  95 !el      integer maxhess,maxlight,maxvib
  96 !el      integer maxgeo,maxcell,maxring
  97 !el      integer maxfix,maxbio
  98 !el      integer maxamino,maxnuc,maxbnd
  99 !el      integer maxang,maxtors,maxpi
 100 !el      integer maxpib,maxpit
 101       integer :: maxatm !=2*nres        !maxres2 maxres2=2*maxres
 102       integer,parameter :: maxval=8
 103       integer,parameter :: maxgrp=1000
 104       integer,parameter :: maxtyp=3000
 105       integer,parameter :: maxclass=500
 106       integer,parameter :: maxkey=10000
 107       integer,parameter :: maxrot=1000
 108       integer,parameter :: maxopt=1000
 109       integer,parameter :: maxhess=1000000
 110       integer :: maxlight       !=8*maxatm
 111       integer,parameter :: maxvib=1000
 112       integer,parameter :: maxgeo=1000
 113       integer,parameter :: maxcell=10000
 114       integer,parameter :: maxring=10000
 115       integer,parameter :: maxfix=10000
 116       integer,parameter :: maxbio=10000
 117       integer,parameter :: maxamino=31
 118       integer,parameter :: maxnuc=12
 119       integer :: maxbnd         !=2*maxatm
 120       integer :: maxang         !=3*maxatm
 121       integer :: maxtors        !=4*maxatm
 122       integer,parameter :: maxpi=100
 123       integer,parameter :: maxpib=2*maxpi
 124       integer,parameter :: maxpit=4*maxpi
 125 !-----------------------------------------------------------------------------
 126 ! Maximum number of seed
 127       integer,parameter :: max_seed=1
 128 !-----------------------------------------------------------------------------
 129       real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
 130 !      common /stochcalc/ stochforcvec
 131 !-----------------------------------------------------------------------------
 132 !      common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
 133       real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
 134       real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
 135       real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
 136 !-----------------------------------------------------------------------------
 137 !      common /syfek/ subroutines: friction_force,setup_fricmat
 138 !el      real(kind=8),dimension(:),allocatable :: gamvec        !(MAXRES6) or (MAXRES2)
 139 !-----------------------------------------------------------------------------
 140 !      common /przechowalnia/ subroutines: friction_force,setup_fricmat
 141       real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
 142 !-----------------------------------------------------------------------------
 143 !      common /przechowalnia/ subroutine: setup_fricmat
 144 !el      real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
 145 !-----------------------------------------------------------------------------
 146 !
 147 !
 148 !-----------------------------------------------------------------------------
 149       contains
 150 !-----------------------------------------------------------------------------
 151 ! brown_step.f
 152 !-----------------------------------------------------------------------------
 153       subroutine brown_step(itime)
 154 !------------------------------------------------
 155 !  Perform a single Euler integration step of Brownian dynamics
 156 !------------------------------------------------
 157 !      implicit real*8 (a-h,o-z)
 158       use comm_gucio
 159       use control, only: tcpu
 160       use control_data
 161       use energy_data
 162 !      use io_conf, only:cartprint
 163 !      include 'DIMENSIONS'
 164 #ifdef MPI
 165       include 'mpif.h'
 166 #endif
 167 !      include 'COMMON.CONTROL'
 168 !      include 'COMMON.VAR'
 169 !      include 'COMMON.MD'
 170 !#ifndef LANG0
 171 !      include 'COMMON.LANGEVIN'
 172 !#else
 173 !      include 'COMMON.LANGEVIN.lang0'
 174 !#endif
 175 !      include 'COMMON.CHAIN'
 176 !      include 'COMMON.DERIV'
 177 !      include 'COMMON.GEO'
 178 !      include 'COMMON.LOCAL'
 179 !      include 'COMMON.INTERACT'
 180 !      include 'COMMON.IOUNITS'
 181 !      include 'COMMON.NAMES'
 182 !      include 'COMMON.TIME1'
 183       real(kind=8),dimension(6*nres) :: zapas   !(MAXRES6) maxres6=6*maxres
 184       integer :: rstcount       !ilen,
 185 !el      external ilen
 186 !el      real(kind=8),dimension(6*nres) :: stochforcvec  !(MAXRES6) maxres6=6*maxres
 187       real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat  !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
 188       real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv       !(maxres2,maxres2) maxres2=2*maxres
 189       real(kind=8),dimension(6*nres,6*nres) :: Pmat     !(maxres6,maxres6) maxres6=6*maxres
 190       real(kind=8),dimension(6*nres) :: Td      !(maxres6) maxres6=6*maxres
 191       real(kind=8),dimension(2*nres) :: ppvec   !(maxres2) maxres2=2*maxres
 192 !el      common /stochcalc/ stochforcvec
 193 !el      real(kind=8),dimension(3) :: cm        !el
 194 !el      common /gucio/ cm
 195       integer :: itime
 196       logical :: lprn = .false.,lprn1 = .false.
 197       integer :: maxiter = 5
 198       real(kind=8) :: difftol = 1.0d-5
 199       real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
 200       integer :: i,j,nbond,k,ind,ind1,iter
 201       integer :: nres2,nres6
 202       logical :: osob
 203       nres2=2*nres
 204       nres6=6*nres
 205
 206       if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6))   !(MAXRES6) maxres6=6*maxres
 207
 208       nbond=nct-nnt
 209       do i=nnt,nct
 210         if (itype(i).ne.10) nbond=nbond+1
 211       enddo
 212 !
 213       if (lprn1) then
 214         write (iout,*) "Generalized inverse of fricmat"
 215         call matout(dimen,dimen,nres6,nres6,fricmat)
 216       endif
 217       do i=1,dimen
 218         do j=1,nbond
 219           Bmat(i,j)=0.0d0
 220         enddo
 221       enddo
 222       ind=3
 223       ind1=0
 224       do i=nnt,nct-1
 225         ind1=ind1+1
 226         do j=1,3
 227           Bmat(ind+j,ind1)=dC_norm(j,i)
 228         enddo
 229         ind=ind+3
 230       enddo
 231       do i=nnt,nct
 232         if (itype(i).ne.10) then
 233           ind1=ind1+1
 234           do j=1,3
 235             Bmat(ind+j,ind1)=dC_norm(j,i+nres)
 236           enddo
 237           ind=ind+3
 238         endif
 239       enddo
 240       if (lprn1) then
 241         write (iout,*) "Matrix Bmat"
 242         call MATOUT(nbond,dimen,nres6,nres6,Bmat)
 243       endif
 244       do i=1,dimen
 245         do j=1,nbond
 246           GBmat(i,j)=0.0d0
 247           do k=1,dimen
 248             GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
 249           enddo
 250         enddo
 251       enddo
 252       if (lprn1) then
 253         write (iout,*) "Matrix GBmat"
 254         call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
 255       endif
 256       do i=1,nbond
 257         do j=1,nbond
 258           Cmat_(i,j)=0.0d0
 259           do k=1,dimen
 260             Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
 261           enddo
 262         enddo
 263       enddo
 264       if (lprn1) then
 265         write (iout,*) "Matrix Cmat"
 266         call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
 267       endif
 268       call matinvert(nbond,nres2,Cmat_,Cinv,osob)
 269       if (lprn1) then
 270         write (iout,*) "Matrix Cinv"
 271         call MATOUT(nbond,nbond,nres2,nres2,Cinv)
 272       endif
 273       do i=1,dimen
 274         do j=1,nbond
 275           Tmat(i,j)=0.0d0
 276           do k=1,nbond
 277             Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
 278           enddo
 279         enddo
 280       enddo
 281       if (lprn1) then
 282         write (iout,*) "Matrix Tmat"
 283         call MATOUT(nbond,dimen,nres6,nres2,Tmat)
 284       endif
 285       do i=1,dimen
 286         do j=1,dimen
 287           if (i.eq.j) then
 288             Pmat(i,j)=1.0d0
 289           else
 290             Pmat(i,j)=0.0d0
 291           endif
 292           do k=1,nbond
 293             Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
 294           enddo
 295         enddo
 296       enddo
 297       if (lprn1) then
 298         write (iout,*) "Matrix Pmat"
 299         call MATOUT(dimen,dimen,nres6,nres6,Pmat)
 300       endif
 301       do i=1,dimen
 302         Td(i)=0.0d0
 303         ind=0
 304         do k=nnt,nct-1
 305           ind=ind+1
 306           Td(i)=Td(i)+vbl*Tmat(i,ind)
 307         enddo
 308         do k=nnt,nct
 309           if (itype(k).ne.10) then
 310             ind=ind+1
 311             Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
 312           endif
 313         enddo
 314       enddo
 315       if (lprn1) then
 316         write (iout,*) "Vector Td"
 317         do i=1,dimen
 318           write (iout,'(i5,f10.5)') i,Td(i)
 319         enddo
 320       endif
 321       call stochastic_force(stochforcvec)
 322       if (lprn) then
 323         write (iout,*) "stochforcvec"
 324         do i=1,dimen
 325           write (iout,*) i,stochforcvec(i)
 326         enddo
 327       endif
 328       do j=1,3
 329         zapas(j)=-gcart(j,0)+stochforcvec(j)
 330         d_t_work(j)=d_t(j,0)
 331         dC_work(j)=dC_old(j,0)
 332       enddo
 333       ind=3
 334       do i=nnt,nct-1
 335         do j=1,3
 336           ind=ind+1
 337           zapas(ind)=-gcart(j,i)+stochforcvec(ind)
 338           dC_work(ind)=dC_old(j,i)
 339         enddo
 340       enddo
 341       do i=nnt,nct
 342         if (itype(i).ne.10) then
 343           do j=1,3
 344             ind=ind+1
 345             zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
 346             dC_work(ind)=dC_old(j,i+nres)
 347           enddo
 348         endif
 349       enddo
 350
 351       if (lprn) then
 352         write (iout,*) "Initial d_t_work"
 353         do i=1,dimen
 354           write (iout,*) i,d_t_work(i)
 355         enddo
 356       endif
 357
 358       do i=1,dimen
 359         d_t_work(i)=0.0d0
 360         do j=1,dimen
 361           d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
 362         enddo
 363       enddo
 364
 365       do i=1,dimen
 366         zapas(i)=Td(i)
 367         do j=1,dimen
 368           zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
 369         enddo
 370       enddo
 371       if (lprn1) then
 372         write (iout,*) "Final d_t_work and zapas"
 373         do i=1,dimen
 374           write (iout,*) i,d_t_work(i),zapas(i)
 375         enddo
 376       endif
 377
 378       do j=1,3
 379         d_t(j,0)=d_t_work(j)
 380         dc(j,0)=zapas(j)
 381         dc_work(j)=dc(j,0)
 382       enddo
 383       ind=3
 384       do i=nnt,nct-1
 385         do j=1,3
 386           d_t(j,i)=d_t_work(i)
 387           dc(j,i)=zapas(ind+j)
 388           dc_work(ind+j)=dc(j,i)
 389         enddo
 390         ind=ind+3
 391       enddo
 392       do i=nnt,nct
 393         do j=1,3
 394           d_t(j,i+nres)=d_t_work(ind+j)
 395           dc(j,i+nres)=zapas(ind+j)
 396           dc_work(ind+j)=dc(j,i+nres)
 397         enddo
 398         ind=ind+3
 399       enddo
 400       if (lprn) then
 401         call chainbuild_cart
 402         write (iout,*) "Before correction for rotational lengthening"
 403         write (iout,*) "New coordinates",&
 404         " and differences between actual and standard bond lengths"
 405         ind=0
 406         do i=nnt,nct-1
 407           ind=ind+1
 408           xx=vbld(i+1)-vbl
 409           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
 410               i,(dC(j,i),j=1,3),xx
 411         enddo
 412         do i=nnt,nct
 413           if (itype(i).ne.10) then
 414             ind=ind+1
 415             xx=vbld(i+nres)-vbldsc0(1,itype(i))
 416             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
 417              i,(dC(j,i+nres),j=1,3),xx
 418           endif
 419         enddo
 420       endif
 421 ! Second correction (rotational lengthening)
 422 !      do iter=1,maxiter
 423       diffmax=0.0d0
 424       ind=0
 425       do i=nnt,nct-1
 426         ind=ind+1
 427         blen2 = scalar(dc(1,i),dc(1,i))
 428         ppvec(ind)=2*vbl**2-blen2
 429         diffbond=dabs(vbl-dsqrt(blen2))
 430         if (diffbond.gt.diffmax) diffmax=diffbond
 431         if (ppvec(ind).gt.0.0d0) then
 432           ppvec(ind)=dsqrt(ppvec(ind))
 433         else
 434           ppvec(ind)=0.0d0
 435         endif
 436         if (lprn) then
 437           write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
 438         endif
 439       enddo
 440       do i=nnt,nct
 441         if (itype(i).ne.10) then
 442           ind=ind+1
 443           blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
 444           ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
 445           diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
 446           if (diffbond.gt.diffmax) diffmax=diffbond
 447           if (ppvec(ind).gt.0.0d0) then
 448             ppvec(ind)=dsqrt(ppvec(ind))
 449           else
 450             ppvec(ind)=0.0d0
 451           endif
 452           if (lprn) then
 453             write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
 454           endif
 455         endif
 456       enddo
 457       if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
 458       if (diffmax.lt.difftol) goto 10
 459       do i=1,dimen
 460         Td(i)=0.0d0
 461         do j=1,nbond
 462           Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
 463         enddo
 464       enddo
 465       do i=1,dimen
 466         zapas(i)=Td(i)
 467         do j=1,dimen
 468           zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
 469         enddo
 470       enddo
 471       do j=1,3
 472         dc(j,0)=zapas(j)
 473         dc_work(j)=zapas(j)
 474       enddo
 475       ind=3
 476       do i=nnt,nct-1
 477         do j=1,3
 478           dc(j,i)=zapas(ind+j)
 479           dc_work(ind+j)=zapas(ind+j)
 480         enddo
 481         ind=ind+3
 482       enddo
 483       do i=nnt,nct
 484         if (itype(i).ne.10) then
 485           do j=1,3
 486             dc(j,i+nres)=zapas(ind+j)
 487             dc_work(ind+j)=zapas(ind+j)
 488           enddo
 489           ind=ind+3
 490         endif
 491       enddo
 492 !   Building the chain from the newly calculated coordinates
 493       call chainbuild_cart
 494       if(ntwe.ne.0) then
 495       if (large.and. mod(itime,ntwe).eq.0) then
 496         write (iout,*) "Cartesian and internal coordinates: step 1"
 497         call cartprint
 498         call intout
 499         write (iout,'(a)') "Potential forces"
 500         do i=0,nres
 501           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
 502           (-gxcart(j,i),j=1,3)
 503         enddo
 504         write (iout,'(a)') "Stochastic forces"
 505         do i=0,nres
 506           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
 507           (stochforc(j,i+nres),j=1,3)
 508         enddo
 509         write (iout,'(a)') "Velocities"
 510         do i=0,nres
 511           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
 512           (d_t(j,i+nres),j=1,3)
 513         enddo
 514       endif
 515       endif
 516       if (lprn) then
 517         write (iout,*) "After correction for rotational lengthening"
 518         write (iout,*) "New coordinates",&
 519         " and differences between actual and standard bond lengths"
 520         ind=0
 521         do i=nnt,nct-1
 522           ind=ind+1
 523           xx=vbld(i+1)-vbl
 524           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
 525               i,(dC(j,i),j=1,3),xx
 526         enddo
 527         do i=nnt,nct
 528           if (itype(i).ne.10) then
 529             ind=ind+1
 530             xx=vbld(i+nres)-vbldsc0(1,itype(i))
 531             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
 532              i,(dC(j,i+nres),j=1,3),xx
 533           endif
 534         enddo
 535       endif
 536 !      ENDDO
 537 !      write (iout,*) "Too many attempts at correcting the bonds"
 538 !      stop
 539    10 continue
 540 #ifdef MPI
 541       tt0 =MPI_Wtime()
 542 #else
 543       tt0 = tcpu()
 544 #endif
 545 ! Calculate energy and forces
 546       call zerograd
 547       call etotal(potEcomp)
 548       potE=potEcomp(0)-potEcomp(20)
 549       call cartgrad
 550       totT=totT+d_time
 551 !  Calculate the kinetic and total energy and the kinetic temperature
 552       call kinetic(EK)
 553 #ifdef MPI
 554       t_enegrad=t_enegrad+MPI_Wtime()-tt0
 555 #else
 556       t_enegrad=t_enegrad+tcpu()-tt0
 557 #endif
 558       totE=EK+potE
 559       kinetic_T=2.0d0/(dimen*Rb)*EK
 560       return
 561       end subroutine brown_step
 562 !-----------------------------------------------------------------------------
 563 ! gauss.f
 564 !-----------------------------------------------------------------------------
 565       subroutine gauss(RO,AP,MT,M,N,*)
 566 !
 567 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
 568 ! RO IS A SQUARE MATRIX
 569 ! THE CALCULATED PRODUCT IS STORED IN AP
 570 ! ABNORMAL EXIT IF RO IS SINGULAR
 571 !
 572       integer :: MT, M, N, M1,I,J,IM,&
 573                  I1,MI,MI1
 574       real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
 575       integer :: k
 576 !      real(kind=8) ::
 577
 578       if(M.ne.1)goto 10
 579       X=RO(1,1)
 580       if(dabs(X).le.1.0D-13) return 1
 581       X=1.0/X
 582       do 16 I=1,N
 583 16     AP(1,I)=AP(1,I)*X
 584        return
 585 10     continue
 586         M1=M-1
 587         DO 1 I=1,M1
 588         IM=I
 589         RM=DABS(RO(I,I))
 590         I1=I+1
 591         do 2 J=I1,M
 592         if(DABS(RO(J,I)).LE.RM) goto 2
 593         RM=DABS(RO(J,I))
 594         IM=J
 595 2       continue
 596         If(IM.eq.I)goto 17
 597         do 3 J=1,N
 598         PR=AP(I,J)
 599         AP(I,J)=AP(IM,J)
 600 3       AP(IM,J)=PR
 601         do 4 J=I,M
 602         PR=RO(I,J)
 603         RO(I,J)=RO(IM,J)
 604 4       RO(IM,J)=PR
 605 17      X=RO(I,I)
 606         if(dabs(X).le.1.0E-13) return 1
 607         X=1.0/X
 608         do 5 J=1,N
 609 5       AP(I,J)=X*AP(I,J)
 610         do 6 J=I1,M
 611 6       RO(I,J)=X*RO(I,J)
 612         do 7 J=I1,M
 613         Y=RO(J,I)
 614         do 8 K=1,N
 615 8       AP(J,K)=AP(J,K)-Y*AP(I,K)
 616         do 9 K=I1,M
 617 9       RO(J,K)=RO(J,K)-Y*RO(I,K)
 618 7       continue
 619 1       continue
 620         X=RO(M,M)
 621         if(dabs(X).le.1.0E-13) return 1
 622         X=1.0/X
 623         do 11 J=1,N
 624 11      AP(M,J)=X*AP(M,J)
 625         do 12 I=1,M1
 626         MI=M-I
 627         MI1=MI+1
 628         do 14 J=1,N
 629         X=AP(MI,J)
 630         do 15 K=MI1,M
 631 15      X=X-AP(K,J)*RO(MI,K)
 632 14      AP(MI,J)=X
 633 12      continue
 634       return
 635       end subroutine gauss
 636 !-----------------------------------------------------------------------------
 637 ! kinetic_lesyng.f
 638 !-----------------------------------------------------------------------------
 639       subroutine kinetic(KE_total)
 640 !----------------------------------------------------------------
 641 !   This subroutine calculates the total kinetic energy of the chain
 642 !-----------------------------------------------------------------
 643       use energy_data
 644 !      implicit real*8 (a-h,o-z)
 645 !      include 'DIMENSIONS'
 646 !      include 'COMMON.VAR'
 647 !      include 'COMMON.CHAIN'
 648 !      include 'COMMON.DERIV'
 649 !      include 'COMMON.GEO'
 650 !      include 'COMMON.LOCAL'
 651 !      include 'COMMON.INTERACT'
 652 !      include 'COMMON.MD'
 653 !      include 'COMMON.IOUNITS'
 654       real(kind=8) :: KE_total
 655
 656       integer :: i,j,k,iti
 657       real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
 658        mag1,mag2,v(3)
 659
 660       KEt_p=0.0d0
 661       KEt_sc=0.0d0
 662 !      write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
 663 !   The translational part for peptide virtual bonds
 664       do j=1,3
 665         incr(j)=d_t(j,0)
 666       enddo
 667       do i=nnt,nct-1
 668 !        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
 669         do j=1,3
 670           v(j)=incr(j)+0.5d0*d_t(j,i)
 671         enddo
 672         vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
 673         KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
 674         do j=1,3
 675           incr(j)=incr(j)+d_t(j,i)
 676         enddo
 677       enddo
 678 !      write(iout,*) 'KEt_p', KEt_p
 679 ! The translational part for the side chain virtual bond
 680 ! Only now we can initialize incr with zeros. It must be equal
 681 ! to the velocities of the first Calpha.
 682       do j=1,3
 683         incr(j)=d_t(j,0)
 684       enddo
 685       do i=nnt,nct
 686         iti=iabs(itype(i))
 687         if (itype(i).eq.10) then
 688           do j=1,3
 689             v(j)=incr(j)
 690           enddo
 691         else
 692           do j=1,3
 693             v(j)=incr(j)+d_t(j,nres+i)
 694           enddo
 695         endif
 696 !        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
 697 !        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
 698         KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
 699         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
 700         do j=1,3
 701           incr(j)=incr(j)+d_t(j,i)
 702         enddo
 703       enddo
 704 !      goto 111
 705 !      write(iout,*) 'KEt_sc', KEt_sc
 706 !  The part due to stretching and rotation of the peptide groups
 707        KEr_p=0.0D0
 708        do i=nnt,nct-1
 709 !        write (iout,*) "i",i
 710 !        write (iout,*) "i",i," mag1",mag1," mag2",mag2
 711         do j=1,3
 712           incr(j)=d_t(j,i)
 713         enddo
 714 !        write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
 715           KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
 716           +incr(3)*incr(3))
 717        enddo
 718 !      goto 111
 719 !       write(iout,*) 'KEr_p', KEr_p
 720 !  The rotational part of the side chain virtual bond
 721        KEr_sc=0.0D0
 722        do i=nnt,nct
 723         iti=iabs(itype(i))
 724         if (itype(i).ne.10) then
 725         do j=1,3
 726           incr(j)=d_t(j,nres+i)
 727         enddo
 728 !        write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
 729         KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
 730           incr(3)*incr(3))
 731         endif
 732        enddo
 733 ! The total kinetic energy
 734   111  continue
 735 !       write(iout,*) 'KEr_sc', KEr_sc
 736        KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
 737 !       write (iout,*) "KE_total",KE_total
 738       return
 739       end subroutine kinetic
 740 !-----------------------------------------------------------------------------
 741 ! MD_A-MTS.F
 742 !-----------------------------------------------------------------------------
 743       subroutine MD
 744 !------------------------------------------------
 745 !  The driver for molecular dynamics subroutines
 746 !------------------------------------------------
 747       use comm_gucio
 748 !     use MPI
 749       use control, only:tcpu,ovrtim
 750 !      use io_comm, only:ilen
 751       use control_data
 752       use compare, only:secondary2,hairpin
 753       use io, only:cartout,statout
 754 !      implicit real*8 (a-h,o-z)
 755 !      include 'DIMENSIONS'
 756 #ifdef MPI
 757       include "mpif.h"
 758       integer :: IERROR,ERRCODE
 759 #endif
 760 !      include 'COMMON.SETUP'
 761 !      include 'COMMON.CONTROL'
 762 !      include 'COMMON.VAR'
 763 !      include 'COMMON.MD'
 764 !#ifndef LANG0
 765 !      include 'COMMON.LANGEVIN'
 766 !#else
 767 !      include 'COMMON.LANGEVIN.lang0'
 768 !#endif
 769 !      include 'COMMON.CHAIN'
 770 !      include 'COMMON.DERIV'
 771 !      include 'COMMON.GEO'
 772 !      include 'COMMON.LOCAL'
 773 !      include 'COMMON.INTERACT'
 774 !      include 'COMMON.IOUNITS'
 775 !      include 'COMMON.NAMES'
 776 !      include 'COMMON.TIME1'
 777 !      include 'COMMON.HAIRPIN'
 778       real(kind=8),dimension(3) :: L,vcm
 779 #ifdef VOUT
 780       real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
 781 #endif
 782       integer :: rstcount       !ilen,
 783 !el      external ilen
 784       character(len=50) :: tytul
 785 !el      common /gucio/ cm
 786       integer :: itime,i,j,nharp
 787       integer,dimension(4,nres/3) :: iharp      !(4,nres/3)(4,maxres/3)
 788 !      logical :: ovrtim
 789       real(kind=8) :: tt0,scalfac
 790       integer :: nres2
 791       nres2=2*nres
 792 !
 793 #ifdef MPI
 794       if (ilen(tmpdir).gt.0) &
 795         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
 796               //liczba(:ilen(liczba))//'.rst')
 797 #else
 798       if (ilen(tmpdir).gt.0) &
 799         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
 800 #endif
 801       t_MDsetup=0.0d0
 802       t_langsetup=0.0d0
 803       t_MD=0.0d0
 804       t_enegrad=0.0d0
 805       t_sdsetup=0.0d0
 806       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
 807 #ifdef MPI
 808       tt0=MPI_Wtime()
 809 #else
 810       tt0 = tcpu()
 811 #endif
 812 ! Determine the inverse of the inertia matrix.
 813       call setup_MD_matrices
 814 ! Initialize MD
 815       call init_MD
 816 #ifdef MPI
 817       t_MDsetup = MPI_Wtime()-tt0
 818 #else
 819       t_MDsetup = tcpu()-tt0
 820 #endif
 821       rstcount=0
 822 !   Entering the MD loop
 823 #ifdef MPI
 824       tt0 = MPI_Wtime()
 825 #else
 826       tt0 = tcpu()
 827 #endif
 828       if (lang.eq.2 .or. lang.eq.3) then
 829 #ifndef   LANG0
 830         call setup_fricmat
 831         if (lang.eq.2) then
 832           call sd_verlet_p_setup
 833         else
 834           call sd_verlet_ciccotti_setup
 835         endif
 836         do i=1,dimen
 837           do j=1,dimen
 838             pfric0_mat(i,j,0)=pfric_mat(i,j)
 839             afric0_mat(i,j,0)=afric_mat(i,j)
 840             vfric0_mat(i,j,0)=vfric_mat(i,j)
 841             prand0_mat(i,j,0)=prand_mat(i,j)
 842             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 843             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 844           enddo
 845         enddo
 846         flag_stoch(0)=.true.
 847         do i=1,maxflag_stoch
 848           flag_stoch(i)=.false.
 849         enddo
 850 #else
 851         write (iout,*) &
 852          "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 853 #ifdef MPI
 854         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 855 #endif
 856         stop
 857 #endif
 858       else if (lang.eq.1 .or. lang.eq.4) then
 859         call setup_fricmat
 860       endif
 861 #ifdef MPI
 862       t_langsetup=MPI_Wtime()-tt0
 863       tt0=MPI_Wtime()
 864 #else
 865       t_langsetup=tcpu()-tt0
 866       tt0=tcpu()
 867 #endif
 868       do itime=1,n_timestep
 869         if (ovrtim()) exit
 870         if (large.and. mod(itime,ntwe).eq.0) &
 871           write (iout,*) "itime",itime
 872         rstcount=rstcount+1
 873         if (lang.gt.0 .and. surfarea .and. &
 874             mod(itime,reset_fricmat).eq.0) then
 875           if (lang.eq.2 .or. lang.eq.3) then
 876 #ifndef LANG0
 877             call setup_fricmat
 878             if (lang.eq.2) then
 879               call sd_verlet_p_setup
 880             else
 881               call sd_verlet_ciccotti_setup
 882             endif
 883             do i=1,dimen
 884               do j=1,dimen
 885                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 886                 afric0_mat(i,j,0)=afric_mat(i,j)
 887                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 888                 prand0_mat(i,j,0)=prand_mat(i,j)
 889                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 890                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 891               enddo
 892             enddo
 893             flag_stoch(0)=.true.
 894             do i=1,maxflag_stoch
 895               flag_stoch(i)=.false.
 896             enddo
 897 #endif
 898           else if (lang.eq.1 .or. lang.eq.4) then
 899             call setup_fricmat
 900           endif
 901           write (iout,'(a,i10)') &
 902             "Friction matrix reset based on surface area, itime",itime
 903         endif
 904         if (reset_vel .and. tbf .and. lang.eq.0 &
 905             .and. mod(itime,count_reset_vel).eq.0) then
 906           call random_vel
 907           write(iout,'(a,f20.2)') &
 908            "Velocities reset to random values, time",totT
 909           do i=0,2*nres
 910             do j=1,3
 911               d_t_old(j,i)=d_t(j,i)
 912             enddo
 913           enddo
 914         endif
 915         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 916           call inertia_tensor
 917           call vcm_vel(vcm)
 918           do j=1,3
 919              d_t(j,0)=d_t(j,0)-vcm(j)
 920           enddo
 921           call kinetic(EK)
 922           kinetic_T=2.0d0/(dimen3*Rb)*EK
 923           scalfac=dsqrt(T_bath/kinetic_T)
 924           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 925           do i=0,2*nres
 926             do j=1,3
 927               d_t_old(j,i)=scalfac*d_t(j,i)
 928             enddo
 929           enddo
 930         endif
 931         if (lang.ne.4) then
 932           if (RESPA) then
 933 ! Time-reversible RESPA algorithm
 934 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 935             call RESPA_step(itime)
 936           else
 937 ! Variable time step algorithm.
 938             call velverlet_step(itime)
 939           endif
 940         else
 941 #ifdef BROWN
 942           call brown_step(itime)
 943 #else
 944           print *,"Brown dynamics not here!"
 945 #ifdef MPI
 946           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 947 #endif
 948           stop
 949 #endif
 950         endif
 951         if (ntwe.ne.0) then
 952          if (mod(itime,ntwe).eq.0) call statout(itime)
 953 #ifdef VOUT
 954         do j=1,3
 955           v_work(j)=d_t(j,0)
 956         enddo
 957         ind=3
 958         do i=nnt,nct-1
 959           do j=1,3
 960             ind=ind+1
 961             v_work(ind)=d_t(j,i)
 962           enddo
 963         enddo
 964         do i=nnt,nct
 965           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 966             do j=1,3
 967               ind=ind+1
 968               v_work(ind)=d_t(j,i+nres)
 969             enddo
 970           endif
 971         enddo
 972
 973         write (66,'(80f10.5)') &
 974           ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 975         do i=1,ind
 976           v_transf(i)=0.0d0
 977           do j=1,ind
 978             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 979           enddo
 980            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 981         enddo
 982         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 983 #endif
 984         endif
 985         if (mod(itime,ntwx).eq.0) then
 986           write (tytul,'("time",f8.2)') totT
 987           if(mdpdb) then
 988              call hairpin(.true.,nharp,iharp)
 989              call secondary2(.true.)
 990              call pdbout(potE,tytul,ipdb)
 991           else
 992              call cartout(totT)
 993           endif
 994         endif
 995         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
 996            open(irest2,file=rest2name,status='unknown')
 997            write(irest2,*) totT,EK,potE,totE,t_bath
 998            do i=1,2*nres
 999             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1000            enddo
1001            do i=1,2*nres
1002             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1003            enddo
1004           close(irest2)
1005           rstcount=0
1006         endif
1007       enddo
1008
1009 #ifdef MPI
1010       t_MD=MPI_Wtime()-tt0
1011 #else
1012       t_MD=tcpu()-tt0
1013 #endif
1014       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1015         '  Timing  ',&
1016        'MD calculations setup:',t_MDsetup,&
1017        'Energy & gradient evaluation:',t_enegrad,&
1018        'Stochastic MD setup:',t_langsetup,&
1019        'Stochastic MD step setup:',t_sdsetup,&
1020        'MD steps:',t_MD
1021       write (iout,'(/28(1h=),a25,27(1h=))') &
1022        '  End of MD calculation  '
1023 #ifdef TIMING_ENE
1024       write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1025         " eshort",t_eshort
1026       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1027        " time_fricmatmult",time_fricmatmult," time_fsample ",&
1028        time_fsample
1029 #endif
1030       return
1031       end subroutine MD
1032 !-----------------------------------------------------------------------------
1033       subroutine velverlet_step(itime)
1034 !-------------------------------------------------------------------------------
1035 !  Perform a single velocity Verlet step; the time step can be rescaled if
1036 !  increments in accelerations exceed the threshold
1037 !-------------------------------------------------------------------------------
1038 !      implicit real*8 (a-h,o-z)
1039 !      include 'DIMENSIONS'
1040       use comm_gucio
1041       use control, only:tcpu
1042       use control_data
1043 #ifdef MPI
1044       include 'mpif.h'
1045       integer :: ierror,ierrcode
1046       real(kind=8) :: errcode
1047 #endif
1048 !      include 'COMMON.SETUP'
1049 !      include 'COMMON.VAR'
1050 !      include 'COMMON.MD'
1051 !#ifndef LANG0
1052 !      include 'COMMON.LANGEVIN'
1053 !#else
1054 !      include 'COMMON.LANGEVIN.lang0'
1055 !#endif
1056 !      include 'COMMON.CHAIN'
1057 !      include 'COMMON.DERIV'
1058 !      include 'COMMON.GEO'
1059 !      include 'COMMON.LOCAL'
1060 !      include 'COMMON.INTERACT'
1061 !      include 'COMMON.IOUNITS'
1062 !      include 'COMMON.NAMES'
1063 !      include 'COMMON.TIME1'
1064 !      include 'COMMON.MUCA'
1065       real(kind=8),dimension(3) :: vcm,incr
1066       real(kind=8),dimension(3) :: L
1067       integer :: count,rstcount !ilen,
1068 !el      external ilen
1069       character(len=50) :: tytul
1070       integer :: maxcount_scale = 20
1071 !el      common /gucio/ cm
1072 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1073 !el      common /stochcalc/ stochforcvec
1074       integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1075       logical :: scale
1076       real(kind=8) :: epdrift,tt0,fac_time
1077 !
1078       if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres))  !(MAXRES6) maxres6=6*maxres
1079
1080       scale=.true.
1081       icount_scale=0
1082       if (lang.eq.1) then
1083         call sddir_precalc
1084       else if (lang.eq.2 .or. lang.eq.3) then
1085 #ifndef LANG0
1086         call stochastic_force(stochforcvec)
1087 #else
1088         write (iout,*) &
1089          "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1090 #ifdef MPI
1091         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1092 #endif
1093         stop
1094 #endif
1095       endif
1096       itime_scal=0
1097       do while (scale)
1098         icount_scale=icount_scale+1
1099         if (icount_scale.gt.maxcount_scale) then
1100           write (iout,*) &
1101             "ERROR: too many attempts at scaling down the time step. ",&
1102             "amax=",amax,"epdrift=",epdrift,&
1103             "damax=",damax,"edriftmax=",edriftmax,&
1104             "d_time=",d_time
1105           call flush(iout)
1106 #ifdef MPI
1107           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1108 #endif
1109           stop
1110         endif
1111 ! First step of the velocity Verlet algorithm
1112         if (lang.eq.2) then
1113 #ifndef LANG0
1114           call sd_verlet1
1115 #endif
1116         else if (lang.eq.3) then
1117 #ifndef LANG0
1118           call sd_verlet1_ciccotti
1119 #endif
1120         else if (lang.eq.1) then
1121           call sddir_verlet1
1122         else
1123           call verlet1
1124         endif
1125 ! Build the chain from the newly calculated coordinates
1126         call chainbuild_cart
1127         if (rattle) call rattle1
1128         if (ntwe.ne.0) then
1129         if (large.and. mod(itime,ntwe).eq.0) then
1130           write (iout,*) "Cartesian and internal coordinates: step 1"
1131           call cartprint
1132           call intout
1133           write (iout,*) "dC"
1134           do i=0,nres
1135             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1136             (dc(j,i+nres),j=1,3)
1137           enddo
1138           write (iout,*) "Accelerations"
1139           do i=0,nres
1140             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1141             (d_a(j,i+nres),j=1,3)
1142           enddo
1143           write (iout,*) "Velocities, step 1"
1144           do i=0,nres
1145             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1146             (d_t(j,i+nres),j=1,3)
1147           enddo
1148         endif
1149         endif
1150 #ifdef MPI
1151         tt0 = MPI_Wtime()
1152 #else
1153         tt0 = tcpu()
1154 #endif
1155 ! Calculate energy and forces
1156         call zerograd
1157         call etotal(potEcomp)
1158         if (large.and. mod(itime,ntwe).eq.0) &
1159           call enerprint(potEcomp)
1160 #ifdef TIMING_ENE
1161 #ifdef MPI
1162         t_etotal=t_etotal+MPI_Wtime()-tt0
1163 #else
1164         t_etotal=t_etotal+tcpu()-tt0
1165 #endif
1166 #endif
1167         potE=potEcomp(0)-potEcomp(20)
1168         call cartgrad
1169 ! Get the new accelerations
1170         call lagrangian
1171 #ifdef MPI
1172         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1173 #else
1174         t_enegrad=t_enegrad+tcpu()-tt0
1175 #endif
1176 ! Determine maximum acceleration and scale down the timestep if needed
1177         call max_accel
1178         amax=amax/(itime_scal+1)**2
1179         call predict_edrift(epdrift)
1180         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1181 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1182           scale=.true.
1183           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1184             /dlog(2.0d0)+1
1185           itime_scal=itime_scal+ifac_time
1186 !          fac_time=dmin1(damax/amax,0.5d0)
1187           fac_time=0.5d0**ifac_time
1188           d_time=d_time*fac_time
1189           if (lang.eq.2 .or. lang.eq.3) then
1190 #ifndef LANG0
1191 !            write (iout,*) "Calling sd_verlet_setup: 1"
1192 ! Rescale the stochastic forces and recalculate or restore
1193 ! the matrices of tinker integrator
1194             if (itime_scal.gt.maxflag_stoch) then
1195               if (large) write (iout,'(a,i5,a)') &
1196                "Calculate matrices for stochastic step;",&
1197                " itime_scal ",itime_scal
1198               if (lang.eq.2) then
1199                 call sd_verlet_p_setup
1200               else
1201                 call sd_verlet_ciccotti_setup
1202               endif
1203               write (iout,'(2a,i3,a,i3,1h.)') &
1204                "Warning: cannot store matrices for stochastic",&
1205                " integration because the index",itime_scal,&
1206                " is greater than",maxflag_stoch
1207               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1208                " integration Langevin algorithm for better efficiency."
1209             else if (flag_stoch(itime_scal)) then
1210               if (large) write (iout,'(a,i5,a,l1)') &
1211                "Restore matrices for stochastic step; itime_scal ",&
1212                itime_scal," flag ",flag_stoch(itime_scal)
1213               do i=1,dimen
1214                 do j=1,dimen
1215                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1216                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1217                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1218                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1219                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1220                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1221                 enddo
1222               enddo
1223             else
1224               if (large) write (iout,'(2a,i5,a,l1)') &
1225                "Calculate & store matrices for stochastic step;",&
1226                " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1227               if (lang.eq.2) then
1228                 call sd_verlet_p_setup
1229               else
1230                 call sd_verlet_ciccotti_setup
1231               endif
1232               flag_stoch(ifac_time)=.true.
1233               do i=1,dimen
1234                 do j=1,dimen
1235                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1236                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1237                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1238                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1239                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1240                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1241                 enddo
1242               enddo
1243             endif
1244             fac_time=1.0d0/dsqrt(fac_time)
1245             do i=1,dimen
1246               stochforcvec(i)=fac_time*stochforcvec(i)
1247             enddo
1248 #endif
1249           else if (lang.eq.1) then
1250 ! Rescale the accelerations due to stochastic forces
1251             fac_time=1.0d0/dsqrt(fac_time)
1252             do i=1,dimen
1253               d_as_work(i)=d_as_work(i)*fac_time
1254             enddo
1255           endif
1256           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1257             "itime",itime," Timestep scaled down to ",&
1258             d_time," ifac_time",ifac_time," itime_scal",itime_scal
1259         else
1260 ! Second step of the velocity Verlet algorithm
1261           if (lang.eq.2) then
1262 #ifndef LANG0
1263             call sd_verlet2
1264 #endif
1265           else if (lang.eq.3) then
1266 #ifndef LANG0
1267             call sd_verlet2_ciccotti
1268 #endif
1269           else if (lang.eq.1) then
1270             call sddir_verlet2
1271           else
1272             call verlet2
1273           endif
1274           if (rattle) call rattle2
1275           totT=totT+d_time
1276           if (d_time.ne.d_time0) then
1277             d_time=d_time0
1278 #ifndef   LANG0
1279             if (lang.eq.2 .or. lang.eq.3) then
1280               if (large) write (iout,'(a)') &
1281                "Restore original matrices for stochastic step"
1282 !              write (iout,*) "Calling sd_verlet_setup: 2"
1283 ! Restore the matrices of tinker integrator if the time step has been restored
1284               do i=1,dimen
1285                 do j=1,dimen
1286                   pfric_mat(i,j)=pfric0_mat(i,j,0)
1287                   afric_mat(i,j)=afric0_mat(i,j,0)
1288                   vfric_mat(i,j)=vfric0_mat(i,j,0)
1289                   prand_mat(i,j)=prand0_mat(i,j,0)
1290                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1291                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1292                 enddo
1293               enddo
1294             endif
1295 #endif
1296           endif
1297           scale=.false.
1298         endif
1299       enddo
1300 ! Calculate the kinetic and the total energy and the kinetic temperature
1301       call kinetic(EK)
1302       totE=EK+potE
1303 ! diagnostics
1304 !      call kinetic1(EK1)
1305 !      write (iout,*) "step",itime," EK",EK," EK1",EK1
1306 ! end diagnostics
1307 ! Couple the system to Berendsen bath if needed
1308       if (tbf .and. lang.eq.0) then
1309         call verlet_bath
1310       endif
1311       kinetic_T=2.0d0/(dimen3*Rb)*EK
1312 ! Backup the coordinates, velocities, and accelerations
1313       do i=0,2*nres
1314         do j=1,3
1315           dc_old(j,i)=dc(j,i)
1316           d_t_old(j,i)=d_t(j,i)
1317           d_a_old(j,i)=d_a(j,i)
1318         enddo
1319       enddo
1320       if (ntwe.ne.0) then
1321       if (mod(itime,ntwe).eq.0 .and. large) then
1322         write (iout,*) "Velocities, step 2"
1323         do i=0,nres
1324           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1325             (d_t(j,i+nres),j=1,3)
1326         enddo
1327       endif
1328       endif
1329       return
1330       end subroutine velverlet_step
1331 !-----------------------------------------------------------------------------
1332       subroutine RESPA_step(itime)
1333 !-------------------------------------------------------------------------------
1334 !  Perform a single RESPA step.
1335 !-------------------------------------------------------------------------------
1336 !      implicit real*8 (a-h,o-z)
1337 !      include 'DIMENSIONS'
1338       use comm_gucio
1339       use comm_cipiszcze
1340 !     use MPI
1341       use control, only:tcpu
1342       use control_data
1343 !      use io_conf, only:cartprint
1344 #ifdef MPI
1345       include 'mpif.h'
1346       integer :: IERROR,ERRCODE
1347 #endif
1348 !      include 'COMMON.SETUP'
1349 !      include 'COMMON.CONTROL'
1350 !      include 'COMMON.VAR'
1351 !      include 'COMMON.MD'
1352 !#ifndef LANG0
1353 !      include 'COMMON.LANGEVIN'
1354 !#else
1355 !      include 'COMMON.LANGEVIN.lang0'
1356 !#endif
1357 !      include 'COMMON.CHAIN'
1358 !      include 'COMMON.DERIV'
1359 !      include 'COMMON.GEO'
1360 !      include 'COMMON.LOCAL'
1361 !      include 'COMMON.INTERACT'
1362 !      include 'COMMON.IOUNITS'
1363 !      include 'COMMON.NAMES'
1364 !      include 'COMMON.TIME1'
1365       real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1366       real(kind=8),dimension(3) :: L,vcm,incr
1367       real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1368       logical :: PRINT_AMTS_MSG = .false.
1369       integer :: count,rstcount !ilen,
1370 !el      external ilen
1371       character(len=50) :: tytul
1372       integer :: maxcount_scale = 10
1373 !el      common /gucio/ cm
1374 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1375 !el      common /stochcalc/ stochforcvec
1376       integer :: itime,itt,i,j,itsplit
1377       logical :: scale
1378 !el      common /cipiszcze/ itt
1379
1380       real(kind=8) :: epdrift,tt0,epdriftmax
1381       itt = itt_comm
1382
1383       if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres))  !(MAXRES6) maxres6=6*maxres
1384
1385       itt=itime
1386       if (ntwe.ne.0) then
1387       if (large.and. mod(itime,ntwe).eq.0) then
1388         write (iout,*) "***************** RESPA itime",itime
1389         write (iout,*) "Cartesian and internal coordinates: step 0"
1390 !        call cartprint
1391         call pdbout(0.0d0,"cipiszcze",iout)
1392         call intout
1393         write (iout,*) "Accelerations from long-range forces"
1394         do i=0,nres
1395           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1396             (d_a(j,i+nres),j=1,3)
1397         enddo
1398         write (iout,*) "Velocities, step 0"
1399         do i=0,nres
1400           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1401             (d_t(j,i+nres),j=1,3)
1402         enddo
1403       endif
1404       endif
1405 !
1406 ! Perform the initial RESPA step (increment velocities)
1407 !      write (iout,*) "*********************** RESPA ini"
1408       call RESPA_vel
1409       if (ntwe.ne.0) then
1410       if (mod(itime,ntwe).eq.0 .and. large) then
1411         write (iout,*) "Velocities, end"
1412         do i=0,nres
1413           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1414             (d_t(j,i+nres),j=1,3)
1415         enddo
1416       endif
1417       endif
1418 ! Compute the short-range forces
1419 #ifdef MPI
1420       tt0 =MPI_Wtime()
1421 #else
1422       tt0 = tcpu()
1423 #endif
1424 ! 7/2/2009 commented out
1425 !      call zerograd
1426 !      call etotal_short(energia_short)
1427 !      call cartgrad
1428 !      call lagrangian
1429 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1430         do i=0,2*nres
1431           do j=1,3
1432             d_a(j,i)=d_a_short(j,i)
1433           enddo
1434         enddo
1435       if (ntwe.ne.0) then
1436       if (large.and. mod(itime,ntwe).eq.0) then
1437         write (iout,*) "energia_short",energia_short(0)
1438         write (iout,*) "Accelerations from short-range forces"
1439         do i=0,nres
1440           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1441             (d_a(j,i+nres),j=1,3)
1442         enddo
1443       endif
1444       endif
1445 #ifdef MPI
1446         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1447 #else
1448         t_enegrad=t_enegrad+tcpu()-tt0
1449 #endif
1450       do i=0,2*nres
1451         do j=1,3
1452           dc_old(j,i)=dc(j,i)
1453           d_t_old(j,i)=d_t(j,i)
1454           d_a_old(j,i)=d_a(j,i)
1455         enddo
1456       enddo
1457 ! 6/30/08 A-MTS: attempt at increasing the split number
1458       do i=0,2*nres
1459         do j=1,3
1460           dc_old0(j,i)=dc_old(j,i)
1461           d_t_old0(j,i)=d_t_old(j,i)
1462           d_a_old0(j,i)=d_a_old(j,i)
1463         enddo
1464       enddo
1465       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1466       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1467 !
1468       scale=.true.
1469       d_time0=d_time
1470       do while (scale)
1471
1472       scale=.false.
1473 !      write (iout,*) "itime",itime," ntime_split",ntime_split
1474 ! Split the time step
1475       d_time=d_time0/ntime_split
1476 ! Perform the short-range RESPA steps (velocity Verlet increments of
1477 ! positions and velocities using short-range forces)
1478 !      write (iout,*) "*********************** RESPA split"
1479       do itsplit=1,ntime_split
1480         if (lang.eq.1) then
1481           call sddir_precalc
1482         else if (lang.eq.2 .or. lang.eq.3) then
1483 #ifndef LANG0
1484           call stochastic_force(stochforcvec)
1485 #else
1486           write (iout,*) &
1487             "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1488 #ifdef MPI
1489           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1490 #endif
1491           stop
1492 #endif
1493         endif
1494 ! First step of the velocity Verlet algorithm
1495         if (lang.eq.2) then
1496 #ifndef LANG0
1497           call sd_verlet1
1498 #endif
1499         else if (lang.eq.3) then
1500 #ifndef LANG0
1501           call sd_verlet1_ciccotti
1502 #endif
1503         else if (lang.eq.1) then
1504           call sddir_verlet1
1505         else
1506           call verlet1
1507         endif
1508 ! Build the chain from the newly calculated coordinates
1509         call chainbuild_cart
1510         if (rattle) call rattle1
1511         if (ntwe.ne.0) then
1512         if (large.and. mod(itime,ntwe).eq.0) then
1513           write (iout,*) "***** ITSPLIT",itsplit
1514           write (iout,*) "Cartesian and internal coordinates: step 1"
1515           call pdbout(0.0d0,"cipiszcze",iout)
1516 !          call cartprint
1517           call intout
1518           write (iout,*) "Velocities, step 1"
1519           do i=0,nres
1520             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1521               (d_t(j,i+nres),j=1,3)
1522           enddo
1523         endif
1524         endif
1525 #ifdef MPI
1526         tt0 = MPI_Wtime()
1527 #else
1528         tt0 = tcpu()
1529 #endif
1530 ! Calculate energy and forces
1531         call zerograd
1532         call etotal_short(energia_short)
1533         if (large.and. mod(itime,ntwe).eq.0) &
1534           call enerprint(energia_short)
1535 #ifdef TIMING_ENE
1536 #ifdef MPI
1537         t_eshort=t_eshort+MPI_Wtime()-tt0
1538 #else
1539         t_eshort=t_eshort+tcpu()-tt0
1540 #endif
1541 #endif
1542         call cartgrad
1543 ! Get the new accelerations
1544         call lagrangian
1545 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1546         do i=0,2*nres
1547           do j=1,3
1548             d_a_short(j,i)=d_a(j,i)
1549           enddo
1550         enddo
1551         if (ntwe.ne.0) then
1552         if (large.and. mod(itime,ntwe).eq.0) then
1553           write (iout,*)"energia_short",energia_short(0)
1554           write (iout,*) "Accelerations from short-range forces"
1555           do i=0,nres
1556             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1557               (d_a(j,i+nres),j=1,3)
1558           enddo
1559         endif
1560         endif
1561 ! 6/30/08 A-MTS
1562 ! Determine maximum acceleration and scale down the timestep if needed
1563         call max_accel
1564         amax=amax/ntime_split**2
1565         call predict_edrift(epdrift)
1566         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1567          write (iout,*) "amax",amax," damax",damax,&
1568          " epdrift",epdrift," epdriftmax",epdriftmax
1569 ! Exit loop and try with increased split number if the change of
1570 ! acceleration is too big
1571         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1572           if (ntime_split.lt.maxtime_split) then
1573             scale=.true.
1574             ntime_split=ntime_split*2
1575             do i=0,2*nres
1576               do j=1,3
1577                 dc_old(j,i)=dc_old0(j,i)
1578                 d_t_old(j,i)=d_t_old0(j,i)
1579                 d_a_old(j,i)=d_a_old0(j,i)
1580               enddo
1581             enddo
1582             if (PRINT_AMTS_MSG) then
1583             write (iout,*) "acceleration/energy drift too large",amax,&
1584             epdrift," split increased to ",ntime_split," itime",itime,&
1585              " itsplit",itsplit
1586             endif
1587             exit
1588           else
1589             write (iout,*) &
1590             "Uh-hu. Bumpy landscape. Maximum splitting number",&
1591              maxtime_split,&
1592             " already reached!!! Trying to carry on!"
1593           endif
1594         endif
1595 #ifdef MPI
1596         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1597 #else
1598         t_enegrad=t_enegrad+tcpu()-tt0
1599 #endif
1600 ! Second step of the velocity Verlet algorithm
1601         if (lang.eq.2) then
1602 #ifndef LANG0
1603           call sd_verlet2
1604 #endif
1605         else if (lang.eq.3) then
1606 #ifndef LANG0
1607           call sd_verlet2_ciccotti
1608 #endif
1609         else if (lang.eq.1) then
1610           call sddir_verlet2
1611         else
1612           call verlet2
1613         endif
1614         if (rattle) call rattle2
1615 ! Backup the coordinates, velocities, and accelerations
1616         do i=0,2*nres
1617           do j=1,3
1618             dc_old(j,i)=dc(j,i)
1619             d_t_old(j,i)=d_t(j,i)
1620             d_a_old(j,i)=d_a(j,i)
1621           enddo
1622         enddo
1623       enddo
1624
1625       enddo ! while scale
1626
1627 ! Restore the time step
1628       d_time=d_time0
1629 ! Compute long-range forces
1630 #ifdef MPI
1631       tt0 =MPI_Wtime()
1632 #else
1633       tt0 = tcpu()
1634 #endif
1635       call zerograd
1636       call etotal_long(energia_long)
1637       if (large.and. mod(itime,ntwe).eq.0) &
1638           call enerprint(energia_long)
1639 #ifdef TIMING_ENE
1640 #ifdef MPI
1641         t_elong=t_elong+MPI_Wtime()-tt0
1642 #else
1643         t_elong=t_elong+tcpu()-tt0
1644 #endif
1645 #endif
1646       call cartgrad
1647       call lagrangian
1648 #ifdef MPI
1649         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1650 #else
1651         t_enegrad=t_enegrad+tcpu()-tt0
1652 #endif
1653 ! Compute accelerations from long-range forces
1654       if (ntwe.ne.0) then
1655       if (large.and. mod(itime,ntwe).eq.0) then
1656         write (iout,*) "energia_long",energia_long(0)
1657         write (iout,*) "Cartesian and internal coordinates: step 2"
1658 !        call cartprint
1659         call pdbout(0.0d0,"cipiszcze",iout)
1660         call intout
1661         write (iout,*) "Accelerations from long-range forces"
1662         do i=0,nres
1663           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1664             (d_a(j,i+nres),j=1,3)
1665         enddo
1666         write (iout,*) "Velocities, step 2"
1667         do i=0,nres
1668           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1669             (d_t(j,i+nres),j=1,3)
1670         enddo
1671       endif
1672       endif
1673 ! Compute the final RESPA step (increment velocities)
1674 !      write (iout,*) "*********************** RESPA fin"
1675       call RESPA_vel
1676 ! Compute the complete potential energy
1677       do i=0,n_ene
1678         potEcomp(i)=energia_short(i)+energia_long(i)
1679       enddo
1680       potE=potEcomp(0)-potEcomp(20)
1681 !      potE=energia_short(0)+energia_long(0)
1682       totT=totT+d_time
1683 ! Calculate the kinetic and the total energy and the kinetic temperature
1684       call kinetic(EK)
1685       totE=EK+potE
1686 ! Couple the system to Berendsen bath if needed
1687       if (tbf .and. lang.eq.0) then
1688         call verlet_bath
1689       endif
1690       kinetic_T=2.0d0/(dimen3*Rb)*EK
1691 ! Backup the coordinates, velocities, and accelerations
1692       if (ntwe.ne.0) then
1693       if (mod(itime,ntwe).eq.0 .and. large) then
1694         write (iout,*) "Velocities, end"
1695         do i=0,nres
1696           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1697             (d_t(j,i+nres),j=1,3)
1698         enddo
1699       endif
1700       endif
1701       return
1702       end subroutine RESPA_step
1703 !-----------------------------------------------------------------------------
1704       subroutine RESPA_vel
1705 !  First and last RESPA step (incrementing velocities using long-range
1706 !  forces).
1707       use energy_data
1708 !      implicit real*8 (a-h,o-z)
1709 !      include 'DIMENSIONS'
1710 !      include 'COMMON.CONTROL'
1711 !      include 'COMMON.VAR'
1712 !      include 'COMMON.MD'
1713 !      include 'COMMON.CHAIN'
1714 !      include 'COMMON.DERIV'
1715 !      include 'COMMON.GEO'
1716 !      include 'COMMON.LOCAL'
1717 !      include 'COMMON.INTERACT'
1718 !      include 'COMMON.IOUNITS'
1719 !      include 'COMMON.NAMES'
1720       integer :: i,j,inres
1721
1722       do j=1,3
1723         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1724       enddo
1725       do i=nnt,nct-1
1726         do j=1,3
1727           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1728         enddo
1729       enddo
1730       do i=nnt,nct
1731         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1732           inres=i+nres
1733           do j=1,3
1734             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1735           enddo
1736         endif
1737       enddo
1738       return
1739       end subroutine RESPA_vel
1740 !-----------------------------------------------------------------------------
1741       subroutine verlet1
1742 ! Applying velocity Verlet algorithm - step 1 to coordinates
1743       use energy_data
1744 !      implicit real*8 (a-h,o-z)
1745 !      include 'DIMENSIONS'
1746 !      include 'COMMON.CONTROL'
1747 !      include 'COMMON.VAR'
1748 !      include 'COMMON.MD'
1749 !      include 'COMMON.CHAIN'
1750 !      include 'COMMON.DERIV'
1751 !      include 'COMMON.GEO'
1752 !      include 'COMMON.LOCAL'
1753 !      include 'COMMON.INTERACT'
1754 !      include 'COMMON.IOUNITS'
1755 !      include 'COMMON.NAMES'
1756       real(kind=8) :: adt,adt2
1757       integer :: i,j,inres
1758
1759 #ifdef DEBUG
1760       write (iout,*) "VELVERLET1 START: DC"
1761       do i=0,nres
1762         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1763          (dc(j,i+nres),j=1,3)
1764       enddo
1765 #endif
1766       do j=1,3
1767         adt=d_a_old(j,0)*d_time
1768         adt2=0.5d0*adt
1769         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1770         d_t_new(j,0)=d_t_old(j,0)+adt2
1771         d_t(j,0)=d_t_old(j,0)+adt
1772       enddo
1773       do i=nnt,nct-1
1774         do j=1,3
1775           adt=d_a_old(j,i)*d_time
1776           adt2=0.5d0*adt
1777           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1778           d_t_new(j,i)=d_t_old(j,i)+adt2
1779           d_t(j,i)=d_t_old(j,i)+adt
1780         enddo
1781       enddo
1782       do i=nnt,nct
1783         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1784           inres=i+nres
1785           do j=1,3
1786             adt=d_a_old(j,inres)*d_time
1787             adt2=0.5d0*adt
1788             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1789             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1790             d_t(j,inres)=d_t_old(j,inres)+adt
1791           enddo
1792         endif
1793       enddo
1794 #ifdef DEBUG
1795       write (iout,*) "VELVERLET1 END: DC"
1796       do i=0,nres
1797         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1798          (dc(j,i+nres),j=1,3)
1799       enddo
1800 #endif
1801       return
1802       end subroutine verlet1
1803 !-----------------------------------------------------------------------------
1804       subroutine verlet2
1805 !  Step 2 of the velocity Verlet algorithm: update velocities
1806       use energy_data
1807 !      implicit real*8 (a-h,o-z)
1808 !      include 'DIMENSIONS'
1809 !      include 'COMMON.CONTROL'
1810 !      include 'COMMON.VAR'
1811 !      include 'COMMON.MD'
1812 !      include 'COMMON.CHAIN'
1813 !      include 'COMMON.DERIV'
1814 !      include 'COMMON.GEO'
1815 !      include 'COMMON.LOCAL'
1816 !      include 'COMMON.INTERACT'
1817 !      include 'COMMON.IOUNITS'
1818 !      include 'COMMON.NAMES'
1819       integer :: i,j,inres
1820
1821       do j=1,3
1822         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1823       enddo
1824       do i=nnt,nct-1
1825         do j=1,3
1826           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1827         enddo
1828       enddo
1829       do i=nnt,nct
1830         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1831           inres=i+nres
1832           do j=1,3
1833             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1834           enddo
1835         endif
1836       enddo
1837       return
1838       end subroutine verlet2
1839 !-----------------------------------------------------------------------------
1840       subroutine sddir_precalc
1841 ! Applying velocity Verlet algorithm - step 1 to coordinates
1842 !      implicit real*8 (a-h,o-z)
1843 !      include 'DIMENSIONS'
1844       use MPI_data
1845       use control_data
1846 #ifdef MPI
1847       include 'mpif.h'
1848 #endif
1849 !      include 'COMMON.CONTROL'
1850 !      include 'COMMON.VAR'
1851 !      include 'COMMON.MD'
1852 !#ifndef LANG0
1853 !      include 'COMMON.LANGEVIN'
1854 !#else
1855 !      include 'COMMON.LANGEVIN.lang0'
1856 !#endif
1857 !      include 'COMMON.CHAIN'
1858 !      include 'COMMON.DERIV'
1859 !      include 'COMMON.GEO'
1860 !      include 'COMMON.LOCAL'
1861 !      include 'COMMON.INTERACT'
1862 !      include 'COMMON.IOUNITS'
1863 !      include 'COMMON.NAMES'
1864 !      include 'COMMON.TIME1'
1865 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1866 !el      common /stochcalc/ stochforcvec
1867       real(kind=8) :: time00
1868 !
1869 ! Compute friction and stochastic forces
1870 !
1871 #ifdef MPI
1872       time00=MPI_Wtime()
1873       call friction_force
1874       time_fric=time_fric+MPI_Wtime()-time00
1875       time00=MPI_Wtime()
1876       call stochastic_force(stochforcvec)
1877       time_stoch=time_stoch+MPI_Wtime()-time00
1878 #endif
1879 !
1880 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1881 ! forces (d_as_work)
1882 !
1883       call ginv_mult(fric_work, d_af_work)
1884       call ginv_mult(stochforcvec, d_as_work)
1885       return
1886       end subroutine sddir_precalc
1887 !-----------------------------------------------------------------------------
1888       subroutine sddir_verlet1
1889 ! Applying velocity Verlet algorithm - step 1 to velocities
1890 !
1891       use energy_data
1892 !      implicit real*8 (a-h,o-z)
1893 !      include 'DIMENSIONS'
1894 !      include 'COMMON.CONTROL'
1895 !      include 'COMMON.VAR'
1896 !      include 'COMMON.MD'
1897 !#ifndef LANG0
1898 !      include 'COMMON.LANGEVIN'
1899 !#else
1900 !      include 'COMMON.LANGEVIN.lang0'
1901 !#endif
1902 !      include 'COMMON.CHAIN'
1903 !      include 'COMMON.DERIV'
1904 !      include 'COMMON.GEO'
1905 !      include 'COMMON.LOCAL'
1906 !      include 'COMMON.INTERACT'
1907 !      include 'COMMON.IOUNITS'
1908 !      include 'COMMON.NAMES'
1909 ! Revised 3/31/05 AL: correlation between random contributions to
1910 ! position and velocity increments included.
1911       real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1912       real(kind=8) :: adt,adt2
1913       integer :: i,j,ind,inres
1914 !
1915 ! Add the contribution from BOTH friction and stochastic force to the
1916 ! coordinates, but ONLY the contribution from the friction forces to velocities
1917 !
1918       do j=1,3
1919         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1920         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1921         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1922         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1923         d_t(j,0)=d_t_old(j,0)+adt
1924       enddo
1925       ind=3
1926       do i=nnt,nct-1
1927         do j=1,3
1928           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1929           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1930           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1931           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1932           d_t(j,i)=d_t_old(j,i)+adt
1933         enddo
1934         ind=ind+3
1935       enddo
1936       do i=nnt,nct
1937         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1938           inres=i+nres
1939           do j=1,3
1940             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1941             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1942             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1943             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1944             d_t(j,inres)=d_t_old(j,inres)+adt
1945           enddo
1946           ind=ind+3
1947         endif
1948       enddo
1949       return
1950       end subroutine sddir_verlet1
1951 !-----------------------------------------------------------------------------
1952       subroutine sddir_verlet2
1953 !  Calculating the adjusted velocities for accelerations
1954 !
1955       use energy_data
1956 !      implicit real*8 (a-h,o-z)
1957 !      include 'DIMENSIONS'
1958 !      include 'COMMON.CONTROL'
1959 !      include 'COMMON.VAR'
1960 !      include 'COMMON.MD'
1961 !#ifndef LANG0
1962 !      include 'COMMON.LANGEVIN'
1963 !#else
1964 !      include 'COMMON.LANGEVIN.lang0'
1965 !#endif
1966 !      include 'COMMON.CHAIN'
1967 !      include 'COMMON.DERIV'
1968 !      include 'COMMON.GEO'
1969 !      include 'COMMON.LOCAL'
1970 !      include 'COMMON.INTERACT'
1971 !      include 'COMMON.IOUNITS'
1972 !      include 'COMMON.NAMES'
1973       real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
1974       real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
1975       integer :: i,j,ind,inres
1976 ! Revised 3/31/05 AL: correlation between random contributions to
1977 ! position and velocity increments included.
1978 ! The correlation coefficients are calculated at low-friction limit.
1979 ! Also, friction forces are now not calculated with new velocities.
1980
1981 !      call friction_force
1982       call stochastic_force(stochforcvec)
1983 !
1984 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1985 ! forces (d_as_work)
1986 !
1987       call ginv_mult(stochforcvec, d_as_work1)
1988
1989 !
1990 ! Update velocities
1991 !
1992       do j=1,3
1993         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
1994           +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1995       enddo
1996       ind=3
1997       do i=nnt,nct-1
1998         do j=1,3
1999           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2000            +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2001         enddo
2002         ind=ind+3
2003       enddo
2004       do i=nnt,nct
2005         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2006           inres=i+nres
2007           do j=1,3
2008             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2009              +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2010              +cos60*d_as_work1(ind+j))*d_time
2011           enddo
2012           ind=ind+3
2013         endif
2014       enddo
2015       return
2016       end subroutine sddir_verlet2
2017 !-----------------------------------------------------------------------------
2018       subroutine max_accel
2019 !
2020 ! Find the maximum difference in the accelerations of the the sites
2021 ! at the beginning and the end of the time step.
2022 !
2023       use energy_data
2024 !      implicit real*8 (a-h,o-z)
2025 !      include 'DIMENSIONS'
2026 !      include 'COMMON.CONTROL'
2027 !      include 'COMMON.VAR'
2028 !      include 'COMMON.MD'
2029 !      include 'COMMON.CHAIN'
2030 !      include 'COMMON.DERIV'
2031 !      include 'COMMON.GEO'
2032 !      include 'COMMON.LOCAL'
2033 !      include 'COMMON.INTERACT'
2034 !      include 'COMMON.IOUNITS'
2035       real(kind=8),dimension(3) :: aux,accel,accel_old
2036       real(kind=8) :: dacc
2037       integer :: i,j
2038
2039       do j=1,3
2040 !        aux(j)=d_a(j,0)-d_a_old(j,0)
2041          accel_old(j)=d_a_old(j,0)
2042          accel(j)=d_a(j,0)
2043       enddo
2044       amax=0.0d0
2045       do i=nnt,nct
2046 ! Backbone
2047         if (i.lt.nct) then
2048 ! 7/3/08 changed to asymmetric difference
2049           do j=1,3
2050 !            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2051             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2052             accel(j)=accel(j)+0.5d0*d_a(j,i)
2053 !            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2054             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2055               dacc=dabs(accel(j)-accel_old(j))
2056 !              write (iout,*) i,dacc
2057               if (dacc.gt.amax) amax=dacc
2058             endif
2059           enddo
2060         endif
2061       enddo
2062 ! Side chains
2063       do j=1,3
2064 !        accel(j)=aux(j)
2065         accel_old(j)=d_a_old(j,0)
2066         accel(j)=d_a(j,0)
2067       enddo
2068       if (nnt.eq.2) then
2069         do j=1,3
2070           accel_old(j)=accel_old(j)+d_a_old(j,1)
2071           accel(j)=accel(j)+d_a(j,1)
2072         enddo
2073       endif
2074       do i=nnt,nct
2075         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2076           do j=1,3
2077 !            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2078             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2079             accel(j)=accel(j)+d_a(j,i+nres)
2080           enddo
2081         endif
2082         do j=1,3
2083 !          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2084           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2085             dacc=dabs(accel(j)-accel_old(j))
2086 !            write (iout,*) "side-chain",i,dacc
2087             if (dacc.gt.amax) amax=dacc
2088           endif
2089         enddo
2090         do j=1,3
2091           accel_old(j)=accel_old(j)+d_a_old(j,i)
2092           accel(j)=accel(j)+d_a(j,i)
2093 !          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2094         enddo
2095       enddo
2096       return
2097       end subroutine max_accel
2098 !-----------------------------------------------------------------------------
2099       subroutine predict_edrift(epdrift)
2100 !
2101 ! Predict the drift of the potential energy
2102 !
2103      use energy_data
2104      use control_data, only: lmuca
2105 !      implicit real*8 (a-h,o-z)
2106 !      include 'DIMENSIONS'
2107 !      include 'COMMON.CONTROL'
2108 !      include 'COMMON.VAR'
2109 !      include 'COMMON.MD'
2110 !      include 'COMMON.CHAIN'
2111 !      include 'COMMON.DERIV'
2112 !      include 'COMMON.GEO'
2113 !      include 'COMMON.LOCAL'
2114 !      include 'COMMON.INTERACT'
2115 !      include 'COMMON.IOUNITS'
2116 !      include 'COMMON.MUCA'
2117       real(kind=8) :: epdrift,epdriftij
2118       integer :: i,j
2119 ! Drift of the potential energy
2120       epdrift=0.0d0
2121       do i=nnt,nct
2122 ! Backbone
2123         if (i.lt.nct) then
2124           do j=1,3
2125             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2126             if (lmuca) epdriftij=epdriftij*factor
2127 !            write (iout,*) "back",i,j,epdriftij
2128             if (epdriftij.gt.epdrift) epdrift=epdriftij
2129           enddo
2130         endif
2131 ! Side chains
2132         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2133           do j=1,3
2134             epdriftij= &
2135              dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2136             if (lmuca) epdriftij=epdriftij*factor
2137 !            write (iout,*) "side",i,j,epdriftij
2138             if (epdriftij.gt.epdrift) epdrift=epdriftij
2139           enddo
2140         endif
2141       enddo
2142       epdrift=0.5d0*epdrift*d_time*d_time
2143 !      write (iout,*) "epdrift",epdrift
2144       return
2145       end subroutine predict_edrift
2146 !-----------------------------------------------------------------------------
2147       subroutine verlet_bath
2148 !
2149 !  Coupling to the thermostat by using the Berendsen algorithm
2150 !
2151       use energy_data
2152 !      implicit real*8 (a-h,o-z)
2153 !      include 'DIMENSIONS'
2154 !      include 'COMMON.CONTROL'
2155 !      include 'COMMON.VAR'
2156 !      include 'COMMON.MD'
2157 !      include 'COMMON.CHAIN'
2158 !      include 'COMMON.DERIV'
2159 !      include 'COMMON.GEO'
2160 !      include 'COMMON.LOCAL'
2161 !      include 'COMMON.INTERACT'
2162 !      include 'COMMON.IOUNITS'
2163 !      include 'COMMON.NAMES'
2164       real(kind=8) :: T_half,fact
2165       integer :: i,j,inres
2166 !
2167       T_half=2.0d0/(dimen3*Rb)*EK
2168       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2169 !      write(iout,*) "T_half", T_half
2170 !      write(iout,*) "EK", EK
2171 !      write(iout,*) "fact", fact
2172       do j=1,3
2173         d_t(j,0)=fact*d_t(j,0)
2174       enddo
2175       do i=nnt,nct-1
2176         do j=1,3
2177           d_t(j,i)=fact*d_t(j,i)
2178         enddo
2179       enddo
2180       do i=nnt,nct
2181         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2182           inres=i+nres
2183           do j=1,3
2184             d_t(j,inres)=fact*d_t(j,inres)
2185           enddo
2186         endif
2187       enddo
2188       return
2189       end subroutine verlet_bath
2190 !-----------------------------------------------------------------------------
2191       subroutine init_MD
2192 !  Set up the initial conditions of a MD simulation
2193       use comm_gucio
2194       use energy_data
2195       use control, only:tcpu
2196 !el      use io_basic, only:ilen
2197       use control_data
2198       use MPI_data
2199       use minimm, only:minim_dc,minimize,sc_move
2200       use io_config, only:readrst
2201       use io, only:statout
2202 !      implicit real*8 (a-h,o-z)
2203 !      include 'DIMENSIONS'
2204 #ifdef MP
2205       include 'mpif.h'
2206       character(len=16) :: form
2207       integer :: IERROR,ERRCODE
2208 #endif
2209 !      include 'COMMON.SETUP'
2210 !      include 'COMMON.CONTROL'
2211 !      include 'COMMON.VAR'
2212 !      include 'COMMON.MD'
2213 !#ifndef LANG0
2214 !      include 'COMMON.LANGEVIN'
2215 !#else
2216 !      include 'COMMON.LANGEVIN.lang0'
2217 !#endif
2218 !      include 'COMMON.CHAIN'
2219 !      include 'COMMON.DERIV'
2220 !      include 'COMMON.GEO'
2221 !      include 'COMMON.LOCAL'
2222 !      include 'COMMON.INTERACT'
2223 !      include 'COMMON.IOUNITS'
2224 !      include 'COMMON.NAMES'
2225 !      include 'COMMON.REMD'
2226       real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2227       real(kind=8),dimension(3) :: vcm,incr,L
2228       real(kind=8) :: xv,sigv,lowb,highb
2229       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
2230       character(len=256) :: qstr
2231 !el      integer ilen
2232 !el      external ilen
2233       character(len=50) :: tytul
2234       logical :: file_exist
2235 !el      common /gucio/ cm
2236       integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
2237       real(kind=8) :: etot,tt0
2238       logical :: fail
2239
2240       d_time0=d_time
2241 !      write(iout,*) "d_time", d_time
2242 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2243 ! if the friction coefficients do not depend on surface area
2244       if (lang.gt.0 .and. .not.surfarea) then
2245         do i=nnt,nct-1
2246           stdforcp(i)=stdfp*dsqrt(gamp)
2247         enddo
2248         do i=nnt,nct
2249           stdforcsc(i)=stdfsc(iabs(itype(i))) &
2250                       *dsqrt(gamsc(iabs(itype(i))))
2251         enddo
2252       endif
2253 ! Open the pdb file for snapshotshots
2254 #ifdef MPI
2255       if(mdpdb) then
2256         if (ilen(tmpdir).gt.0) &
2257           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2258             liczba(:ilen(liczba))//".pdb")
2259         open(ipdb,&
2260         file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2261         //".pdb")
2262       else
2263 #ifdef NOXDR
2264         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2265           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2266             liczba(:ilen(liczba))//".x")
2267         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2268         //".x"
2269 #else
2270         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2271           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2272             liczba(:ilen(liczba))//".cx")
2273         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2274         //".cx"
2275 #endif
2276       endif
2277 #else
2278       if(mdpdb) then
2279          if (ilen(tmpdir).gt.0) &
2280            call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2281          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2282       else
2283          if (ilen(tmpdir).gt.0) &
2284            call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2285          cartname=prefix(:ilen(prefix))//"_MD.cx"
2286       endif
2287 #endif
2288       if (usampl) then
2289         write (qstr,'(256(1h ))')
2290         ipos=1
2291         do i=1,nfrag
2292           iq = qinfrag(i,iset)*10
2293           iw = wfrag(i,iset)/100
2294           if (iw.gt.0) then
2295             if(me.eq.king.or..not.out1file) &
2296              write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2297             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2298             ipos=ipos+7
2299           endif
2300         enddo
2301         do i=1,npair
2302           iq = qinpair(i,iset)*10
2303           iw = wpair(i,iset)/100
2304           if (iw.gt.0) then
2305             if(me.eq.king.or..not.out1file) &
2306              write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2307             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2308             ipos=ipos+7
2309           endif
2310         enddo
2311 !        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2312 #ifdef NOXDR
2313 !        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2314 #else
2315 !        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2316 #endif
2317 !        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2318       endif
2319       icg=1
2320       if (rest) then
2321        if (restart1file) then
2322          if (me.eq.king) &
2323            inquire(file=mremd_rst_name,exist=file_exist)
2324            write (*,*) me," Before broadcast: file_exist",file_exist
2325 #ifdef MPI !el
2326          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2327                 IERR)
2328 #endif !el
2329          write (*,*) me," After broadcast: file_exist",file_exist
2330 !        inquire(file=mremd_rst_name,exist=file_exist)
2331         if(me.eq.king.or..not.out1file) &
2332          write(iout,*) "Initial state read by master and distributed"
2333        else
2334          if (ilen(tmpdir).gt.0) &
2335            call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2336             //liczba(:ilen(liczba))//'.rst')
2337         inquire(file=rest2name,exist=file_exist)
2338        endif
2339        if(file_exist) then
2340          if(.not.restart1file) then
2341            if(me.eq.king.or..not.out1file) &
2342             write(iout,*) "Initial state will be read from file ",&
2343             rest2name(:ilen(rest2name))
2344            call readrst
2345          endif
2346          call rescale_weights(t_bath)
2347        else
2348         if(me.eq.king.or..not.out1file)then
2349          if (restart1file) then
2350           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2351              " does not exist"
2352          else
2353           write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2354              " does not exist"
2355          endif
2356          write(iout,*) "Initial velocities randomly generated"
2357         endif
2358         call random_vel
2359         totT=0.0d0
2360        endif
2361       else
2362 ! Generate initial velocities
2363         if(me.eq.king.or..not.out1file) &
2364          write(iout,*) "Initial velocities randomly generated"
2365         call random_vel
2366         totT=0.0d0
2367       endif
2368 !      rest2name = prefix(:ilen(prefix))//'.rst'
2369       if(me.eq.king.or..not.out1file)then
2370        write (iout,*) "Initial velocities"
2371        do i=0,nres
2372          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2373          (d_t(j,i+nres),j=1,3)
2374        enddo
2375 !  Zeroing the total angular momentum of the system
2376        write(iout,*) "Calling the zero-angular momentum subroutine"
2377       endif
2378       call inertia_tensor
2379 !  Getting the potential energy and forces and velocities and accelerations
2380       call vcm_vel(vcm)
2381 !      write (iout,*) "velocity of the center of the mass:"
2382 !      write (iout,*) (vcm(j),j=1,3)
2383       do j=1,3
2384         d_t(j,0)=d_t(j,0)-vcm(j)
2385       enddo
2386 ! Removing the velocity of the center of mass
2387       call vcm_vel(vcm)
2388       if(me.eq.king.or..not.out1file)then
2389        write (iout,*) "vcm right after adjustment:"
2390        write (iout,*) (vcm(j),j=1,3)
2391       endif
2392       if (.not.rest) then
2393          call chainbuild
2394          if(iranconf.ne.0) then
2395           if (overlapsc) then
2396            print *, 'Calling OVERLAP_SC'
2397            call overlap_sc(fail)
2398           endif
2399           if (searchsc) then
2400            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2401            print *,'SC_move',nft_sc,etot
2402            if(me.eq.king.or..not.out1file) &
2403             write(iout,*) 'SC_move',nft_sc,etot
2404           endif
2405
2406           if(dccart)then
2407            print *, 'Calling MINIM_DC'
2408            call minim_dc(etot,iretcode,nfun)
2409           else
2410            call geom_to_var(nvar,varia)
2411            print *,'Calling MINIMIZE.'
2412            call minimize(etot,varia,iretcode,nfun)
2413            call var_to_geom(nvar,varia)
2414           endif
2415           if(me.eq.king.or..not.out1file) &
2416              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2417          endif
2418       endif
2419       call chainbuild_cart
2420       call kinetic(EK)
2421       if (tbf) then
2422         call verlet_bath
2423       endif
2424       kinetic_T=2.0d0/(dimen3*Rb)*EK
2425       if(me.eq.king.or..not.out1file)then
2426        call cartprint
2427        call intout
2428       endif
2429 #ifdef MPI
2430       tt0=MPI_Wtime()
2431 #else
2432       tt0=tcpu()
2433 #endif
2434       call zerograd
2435       call etotal(potEcomp)
2436       if (large) call enerprint(potEcomp)
2437 #ifdef TIMING_ENE
2438 #ifdef MPI
2439       t_etotal=t_etotal+MPI_Wtime()-tt0
2440 #else
2441       t_etotal=t_etotal+tcpu()-tt0
2442 #endif
2443 #endif
2444       potE=potEcomp(0)
2445       call cartgrad
2446       call lagrangian
2447       call max_accel
2448       if (amax*d_time .gt. dvmax) then
2449         d_time=d_time*dvmax/amax
2450         if(me.eq.king.or..not.out1file) write (iout,*) &
2451          "Time step reduced to",d_time,&
2452          " because of too large initial acceleration."
2453       endif
2454       if(me.eq.king.or..not.out1file)then
2455        write(iout,*) "Potential energy and its components"
2456        call enerprint(potEcomp)
2457 !       write(iout,*) (potEcomp(i),i=0,n_ene)
2458       endif
2459       potE=potEcomp(0)-potEcomp(20)
2460       totE=EK+potE
2461       itime=0
2462       if (ntwe.ne.0) call statout(itime)
2463       if(me.eq.king.or..not.out1file) &
2464         write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2465          " Kinetic energy",EK," Potential energy",potE, &
2466          " Total energy",totE," Maximum acceleration ", &
2467          amax
2468       if (large) then
2469         write (iout,*) "Initial coordinates"
2470         do i=1,nres
2471           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2472           (c(j,i+nres),j=1,3)
2473         enddo
2474         write (iout,*) "Initial dC"
2475         do i=0,nres
2476           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2477           (dc(j,i+nres),j=1,3)
2478         enddo
2479         write (iout,*) "Initial velocities"
2480         write (iout,"(13x,' backbone ',23x,' side chain')")
2481         do i=0,nres
2482           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2483           (d_t(j,i+nres),j=1,3)
2484         enddo
2485         write (iout,*) "Initial accelerations"
2486         do i=0,nres
2487 !          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2488           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2489           (d_a(j,i+nres),j=1,3)
2490         enddo
2491       endif
2492       do i=0,2*nres
2493         do j=1,3
2494           dc_old(j,i)=dc(j,i)
2495           d_t_old(j,i)=d_t(j,i)
2496           d_a_old(j,i)=d_a(j,i)
2497         enddo
2498 !        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2499       enddo
2500       if (RESPA) then
2501 #ifdef MPI
2502         tt0 =MPI_Wtime()
2503 #else
2504         tt0 = tcpu()
2505 #endif
2506         call zerograd
2507         call etotal_short(energia_short)
2508         if (large) call enerprint(potEcomp)
2509 #ifdef TIMING_ENE
2510 #ifdef MPI
2511         t_eshort=t_eshort+MPI_Wtime()-tt0
2512 #else
2513         t_eshort=t_eshort+tcpu()-tt0
2514 #endif
2515 #endif
2516         call cartgrad
2517         call lagrangian
2518         if(.not.out1file .and. large) then
2519           write (iout,*) "energia_long",energia_long(0),&
2520            " energia_short",energia_short(0),&
2521            " total",energia_long(0)+energia_short(0)
2522           write (iout,*) "Initial fast-force accelerations"
2523           do i=0,nres
2524             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2525             (d_a(j,i+nres),j=1,3)
2526           enddo
2527         endif
2528 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2529         do i=0,2*nres
2530           do j=1,3
2531             d_a_short(j,i)=d_a(j,i)
2532           enddo
2533         enddo
2534 #ifdef MPI
2535         tt0=MPI_Wtime()
2536 #else
2537         tt0=tcpu()
2538 #endif
2539         call zerograd
2540         call etotal_long(energia_long)
2541         if (large) call enerprint(potEcomp)
2542 #ifdef TIMING_ENE
2543 #ifdef MPI
2544         t_elong=t_elong+MPI_Wtime()-tt0
2545 #else
2546         t_elong=t_elong+tcpu()-tt0
2547 #endif
2548 #endif
2549         call cartgrad
2550         call lagrangian
2551         if(.not.out1file .and. large) then
2552           write (iout,*) "energia_long",energia_long(0)
2553           write (iout,*) "Initial slow-force accelerations"
2554           do i=0,nres
2555             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2556             (d_a(j,i+nres),j=1,3)
2557           enddo
2558         endif
2559 #ifdef MPI
2560         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2561 #else
2562         t_enegrad=t_enegrad+tcpu()-tt0
2563 #endif
2564       endif
2565       return
2566       end subroutine init_MD
2567 !-----------------------------------------------------------------------------
2568       subroutine random_vel
2569
2570 !      implicit real*8 (a-h,o-z)
2571       use energy_data
2572       use random, only:anorm_distr
2573       use MD_data
2574 !      include 'DIMENSIONS'
2575 !      include 'COMMON.CONTROL'
2576 !      include 'COMMON.VAR'
2577 !      include 'COMMON.MD'
2578 !#ifndef LANG0
2579 !      include 'COMMON.LANGEVIN'
2580 !#else
2581 !      include 'COMMON.LANGEVIN.lang0'
2582 !#endif
2583 !      include 'COMMON.CHAIN'
2584 !      include 'COMMON.DERIV'
2585 !      include 'COMMON.GEO'
2586 !      include 'COMMON.LOCAL'
2587 !      include 'COMMON.INTERACT'
2588 !      include 'COMMON.IOUNITS'
2589 !      include 'COMMON.NAMES'
2590 !      include 'COMMON.TIME1'
2591       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
2592 #ifdef FIVEDIAG
2593        real(kind=8) ,allocatable, dimension(:)  :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2594        real(kind=8) :: sumx
2595 #ifdef DEBUG
2596        real(kind=8) ,allocatable, dimension(:)  :: rsold
2597        real (kind=8),allocatable,dimension(:,:) :: matold
2598 #endif
2599 #endif
2600       integer :: i,j,ii,k,ind,mark,imark
2601 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2602 ! First generate velocities in the eigenspace of the G matrix
2603 !      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2604 !      call flush(iout)
2605       xv=0.0d0
2606       ii=0
2607       do i=1,dimen
2608         do k=1,3
2609           ii=ii+1
2610           sigv=dsqrt((Rb*t_bath)/geigen(i))
2611           lowb=-5*sigv
2612           highb=5*sigv
2613           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2614 #ifdef DEBUG
2615           write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2616             " d_t_work_new",d_t_work_new(ii)
2617 #endif
2618         enddo
2619       enddo
2620 #ifdef DEBUG
2621 ! diagnostics
2622       Ek1=0.0d0
2623       ii=0
2624       do i=1,dimen
2625         do k=1,3
2626           ii=ii+1
2627           Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2628         enddo
2629       enddo
2630       write (iout,*) "Ek from eigenvectors",Ek1
2631       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2632 ! end diagnostics
2633 #endif
2634 #ifdef FIVEDIAG
2635 ! Transform velocities to UNRES coordinate space
2636        allocate (DL1(2*nres))
2637        allocate (DDU1(2*nres))
2638        allocate (DL2(2*nres))
2639        allocate (DDU2(2*nres))
2640        allocate (xsolv(2*nres))
2641        allocate (DML(2*nres))
2642        allocate (rs(2*nres))
2643 #ifdef DEBUG
2644        allocate (rsold(2*nres))
2645        allocate (matold(2*nres,2*nres))
2646        matold=0
2647        matold(1,1)=DMorig(1)
2648        matold(1,2)=DU1orig(1)
2649        matold(1,3)=DU2orig(1)
2650        write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2651
2652       do i=2,dimen-1
2653         do j=1,dimen
2654           if (i.eq.j) then
2655              matold(i,j)=DMorig(i)
2656              matold(i,j-1)=DU1orig(i-1)
2657            if (j.ge.3) then
2658                  matold(i,j-2)=DU2orig(i-2)
2659                endif
2660
2661             endif
2662
2663           enddo
2664           do j=1,dimen-1
2665             if (i.eq.j) then
2666               matold(i,j+1)=DU1orig(i)
2667
2668              end if
2669           enddo
2670           do j=1,dimen-2
2671             if(i.eq.j) then
2672                matold(i,j+2)=DU2orig(i)
2673             endif
2674           enddo
2675        enddo
2676        matold(dimen,dimen)=DMorig(dimen)
2677        matold(dimen,dimen-1)=DU1orig(dimen-1)
2678        matold(dimen,dimen-2)=DU2orig(dimen-2)
2679        write(iout,*) "old gmatrix"
2680        call matout(dimen,dimen,2*nres,2*nres,matold)
2681 #endif
2682       d_t_work=0.0d0
2683       do i=1,dimen
2684 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2685        do imark=dimen,1,-1
2686
2687          do j=1,imark-1
2688            DML(j)=DMorig(j)-geigen(i)
2689          enddo
2690          do j=imark+1,dimen
2691            DML(j-1)=DMorig(j)-geigen(i)
2692          enddo
2693          do j=1,imark-2
2694            DDU1(j)=DU1orig(j)
2695          enddo
2696          DDU1(imark-1)=DU2orig(imark-1)
2697          do j=imark+1,dimen-1
2698            DDU1(j-1)=DU1orig(j)
2699          enddo
2700          do j=1,imark-3
2701            DDU2(j)=DU2orig(j)
2702          enddo
2703          DDU2(imark-2)=0.0d0
2704          DDU2(imark-1)=0.0d0
2705          do j=imark,dimen-3
2706            DDU2(j)=DU2orig(j+1)
2707          enddo
2708          do j=1,dimen-3
2709            DL2(j+2)=DDU2(j)
2710          enddo
2711          do j=1,dimen-2
2712            DL1(j+1)=DDU1(j)
2713          enddo
2714 #ifdef DEBUG
2715          write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2716          write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2717          write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2718          write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2719          write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2720          write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2721 #endif
2722          rs=0.0d0
2723          if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2724          if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2725          if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2726          if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2727 #ifdef DEBUG
2728          rsold=rs
2729 #endif
2730 !         write (iout,*) "Vector rs"
2731 !         do j=1,dimen-1
2732 !           write (iout,*) j,rs(j)
2733 !         enddo
2734          xsolv=0.0d0
2735          call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2736
2737          if (mark.eq.1) then
2738
2739 #ifdef DEBUG
2740 ! Check solution
2741          do j=1,imark-1
2742            sumx=-geigen(i)*xsolv(j)
2743            do k=1,imark-1
2744              sumx=sumx+matold(j,k)*xsolv(k)
2745            enddo
2746            do k=imark+1,dimen
2747              sumx=sumx+matold(j,k)*xsolv(k-1)
2748            enddo
2749            write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2750          enddo
2751          do j=imark+1,dimen
2752            sumx=-geigen(i)*xsolv(j-1)
2753            do k=1,imark-1
2754              sumx=sumx+matold(j,k)*xsolv(k)
2755            enddo
2756            do k=imark+1,dimen
2757              sumx=sumx+matold(j,k)*xsolv(k-1)
2758            enddo
2759            write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2760          enddo
2761 ! end check
2762          write (iout,*)&
2763           "Solution of equations system",i," for eigenvalue",geigen(i)
2764          do j=1,dimen-1
2765            write(iout,'(i5,f10.5)') j,xsolv(j)
2766          enddo
2767 #endif
2768          do j=dimen-1,imark,-1
2769            xsolv(j+1)=xsolv(j)
2770          enddo
2771          xsolv(imark)=1.0d0
2772 #ifdef DEBUG
2773          write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2774          do j=1,dimen
2775            write(iout,'(i5,f10.5)') j,xsolv(j)
2776          enddo
2777 #endif
2778 ! Normalize ith eigenvector
2779          sumx=0.0d0
2780          do j=1,dimen
2781            sumx=sumx+xsolv(j)**2
2782          enddo
2783          sumx=dsqrt(sumx)
2784          do j=1,dimen
2785            xsolv(j)=xsolv(j)/sumx
2786          enddo
2787 #ifdef DEBUG
2788          write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2789          do j=1,dimen
2790            write(iout,'(i5,3f10.5)') j,xsolv(j)
2791          enddo
2792 #endif
2793 ! All done at this point for eigenvector i, exit loop
2794          exit
2795
2796          endif ! mark.eq.1
2797
2798        enddo ! imark
2799
2800        if (mark.ne.1) then
2801          write (iout,*) "Unable to find eigenvector",i
2802        endif
2803
2804 !       write (iout,*) "i=",i
2805        do k=1,3
2806 !         write (iout,*) "k=",k
2807          do j=1,dimen
2808            ind=(j-1)*3+k
2809 !           write(iout,*) "ind",ind," ind1",3*(i-1)+k
2810            d_t_work(ind)=d_t_work(ind) &
2811                             +xsolv(j)*d_t_work_new(3*(i-1)+k)
2812          enddo
2813        enddo
2814       enddo ! i (loop over the eigenvectors)
2815
2816 #ifdef DEBUG
2817       write (iout,*) "d_t_work"
2818       do i=1,3*dimen
2819         write (iout,"(i5,f10.5)") i,d_t_work(i)
2820       enddo
2821       Ek1=0.0d0
2822       ii=0
2823       do i=nnt,nct
2824         if (itype(i).eq.10) then
2825           j=ii+3
2826         else
2827           j=ii+6
2828         endif
2829         if (i.lt.nct) then
2830           do k=1,3
2831 !            write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2832             Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2833             0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2834           enddo
2835         endif
2836         if (itype(i).ne.10) ii=ii+3
2837         write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
2838         write (iout,*) "ii",ii
2839         do k=1,3
2840           ii=ii+1
2841           write (iout,*) "k",k," ii",ii,"EK1",EK1
2842           if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
2843           Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2
2844         enddo
2845         write (iout,*) "i",i," ii",ii
2846       enddo
2847       write (iout,*) "Ek from d_t_work",Ek1
2848       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2849 #endif
2850       deallocate(DDU1)
2851       deallocate(DDU2)
2852       deallocate(DL2)
2853       deallocate(DL1)
2854       deallocate(xsolv)
2855       deallocate(DML)
2856       deallocate(rs)
2857 #ifdef DEBUG
2858       deallocate(matold)
2859       deallocate(rsold)
2860 #endif
2861       ind=1
2862       do j=nnt,nct
2863         do k=1,3
2864           d_t(k,j)=d_t_work(ind)
2865           ind=ind+1
2866         enddo
2867         if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then
2868           do k=1,3
2869             d_t(k,j+nres)=d_t_work(ind)
2870             ind=ind+1
2871           enddo
2872         end if
2873       enddo
2874 #ifdef DEBUG
2875       write (iout,*) "Random velocities in the Calpha,SC space"
2876       do i=nnt,nct
2877         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2878       enddo
2879       do i=nnt,nct
2880         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2881       enddo
2882 #endif
2883       do j=1,3
2884         d_t(j,0)=d_t(j,nnt)
2885       enddo
2886       do i=nnt,nct
2887         if (itype(i).eq.10) then
2888           do j=1,3
2889             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2890           enddo
2891         else
2892           do j=1,3
2893             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2894             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2895           enddo
2896         end if
2897       enddo
2898 #ifdef DEBUG
2899       write (iout,*)"Random velocities in the virtual-bond-vector space"
2900       do i=nnt,nct-1
2901         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2902       enddo
2903       do i=nnt,nct
2904         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2905       enddo
2906       call kinetic(Ek1)
2907       write (iout,*) "Ek from d_t_work",Ek1
2908       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2909 #endif
2910 #else
2911       do k=0,2
2912         do i=1,dimen
2913           ind=(i-1)*3+k+1
2914           d_t_work(ind)=0.0d0
2915           do j=1,dimen
2916             d_t_work(ind)=d_t_work(ind) &
2917                             +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2918           enddo
2919 !          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2920 !          call flush(iout)
2921         enddo
2922       enddo
2923 ! Transfer to the d_t vector
2924       do j=1,3
2925         d_t(j,0)=d_t_work(j)
2926       enddo
2927       ind=3
2928       do i=nnt,nct-1
2929         do j=1,3
2930           ind=ind+1
2931           d_t(j,i)=d_t_work(ind)
2932         enddo
2933       enddo
2934       do i=nnt,nct
2935         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2936           do j=1,3
2937             ind=ind+1
2938             d_t(j,i+nres)=d_t_work(ind)
2939           enddo
2940         endif
2941       enddo
2942 #endif
2943 !      call kinetic(EK)
2944 !      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
2945 !        2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2946 !      call flush(iout)
2947       return
2948 #undef DEBUG
2949       end subroutine random_vel
2950 !-----------------------------------------------------------------------------
2951 #ifndef LANG0
2952       subroutine sd_verlet_p_setup
2953 ! Sets up the parameters of stochastic Verlet algorithm
2954 !      implicit real*8 (a-h,o-z)
2955 !      include 'DIMENSIONS'
2956       use control, only: tcpu
2957       use control_data
2958 #ifdef MPI
2959       include 'mpif.h'
2960 #endif
2961 !      include 'COMMON.CONTROL'
2962 !      include 'COMMON.VAR'
2963 !      include 'COMMON.MD'
2964 !#ifndef LANG0
2965 !      include 'COMMON.LANGEVIN'
2966 !#else
2967 !      include 'COMMON.LANGEVIN.lang0'
2968 !#endif
2969 !      include 'COMMON.CHAIN'
2970 !      include 'COMMON.DERIV'
2971 !      include 'COMMON.GEO'
2972 !      include 'COMMON.LOCAL'
2973 !      include 'COMMON.INTERACT'
2974 !      include 'COMMON.IOUNITS'
2975 !      include 'COMMON.NAMES'
2976 !      include 'COMMON.TIME1'
2977       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
2978       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
2979       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
2980        prand_vec,vrand_vec1,vrand_vec2  !(MAXRES6) maxres6=6*maxres
2981       logical :: lprn = .false.
2982       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
2983       real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
2984                  gdt9,psig,tt0
2985       integer :: i,maxres2
2986 #ifdef MPI
2987       tt0 = MPI_Wtime()
2988 #else
2989       tt0 = tcpu()
2990 #endif
2991 !
2992 ! AL 8/17/04 Code adapted from tinker
2993 !
2994 ! Get the frictional and random terms for stochastic dynamics in the
2995 ! eigenspace of mass-scaled UNRES friction matrix
2996 !
2997       maxres2=2*nres
2998       do i = 1, dimen
2999             gdt = fricgam(i) * d_time
3000 !
3001 ! Stochastic dynamics reduces to simple MD for zero friction
3002 !
3003             if (gdt .le. zero) then
3004                pfric_vec(i) = 1.0d0
3005                vfric_vec(i) = d_time
3006                afric_vec(i) = 0.5d0 * d_time * d_time
3007                prand_vec(i) = 0.0d0
3008                vrand_vec1(i) = 0.0d0
3009                vrand_vec2(i) = 0.0d0
3010 !
3011 ! Analytical expressions when friction coefficient is large
3012 !
3013             else
3014                if (gdt .ge. gdt_radius) then
3015                   egdt = dexp(-gdt)
3016                   pfric_vec(i) = egdt
3017                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3018                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3019                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3020                   vterm = 1.0d0 - egdt**2
3021                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3022 !
3023 ! Use series expansions when friction coefficient is small
3024 !
3025                else
3026                   gdt2 = gdt * gdt
3027                   gdt3 = gdt * gdt2
3028                   gdt4 = gdt2 * gdt2
3029                   gdt5 = gdt2 * gdt3
3030                   gdt6 = gdt3 * gdt3
3031                   gdt7 = gdt3 * gdt4
3032                   gdt8 = gdt4 * gdt4
3033                   gdt9 = gdt4 * gdt5
3034                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3035                                 - gdt5/120.0d0 + gdt6/720.0d0 &
3036                                 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3037                                 - gdt9/362880.0d0) / fricgam(i)**2
3038                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3039                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3040                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3041                              + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3042                              + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3043                              + 127.0d0*gdt9/90720.0d0
3044                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3045                              - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3046                              - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3047                              - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3048                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3049                              - 17.0d0*gdt2/1280.0d0 &
3050                              + 17.0d0*gdt3/6144.0d0 &
3051                              + 40967.0d0*gdt4/34406400.0d0 &
3052                              - 57203.0d0*gdt5/275251200.0d0 &
3053                              - 1429487.0d0*gdt6/13212057600.0d0)
3054                end if
3055 !
3056 ! Compute the scaling factors of random terms for the nonzero friction case
3057 !
3058                ktm = 0.5d0*d_time/fricgam(i)
3059                psig = dsqrt(ktm*pterm) / fricgam(i)
3060                vsig = dsqrt(ktm*vterm)
3061                rhoc = dsqrt(1.0d0 - rho*rho)
3062                prand_vec(i) = psig
3063                vrand_vec1(i) = vsig * rho
3064                vrand_vec2(i) = vsig * rhoc
3065             end if
3066       end do
3067       if (lprn) then
3068       write (iout,*) &
3069         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3070         " vrand_vec2"
3071       do i=1,dimen
3072         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3073             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3074       enddo
3075       endif
3076 !
3077 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3078 !
3079 #ifndef   LANG0
3080       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3081       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3082       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3083       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3084       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3085       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3086 #endif
3087 #ifdef MPI
3088       t_sdsetup=t_sdsetup+MPI_Wtime()
3089 #else
3090       t_sdsetup=t_sdsetup+tcpu()-tt0
3091 #endif
3092       return
3093       end subroutine sd_verlet_p_setup
3094 !-----------------------------------------------------------------------------
3095       subroutine eigtransf1(n,ndim,ab,d,c)
3096
3097 !el      implicit none
3098       integer :: n,ndim
3099       real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3100       integer :: i,j,k
3101       do i=1,n
3102         do j=1,n
3103           c(i,j)=0.0d0
3104           do k=1,n
3105             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3106           enddo
3107         enddo
3108       enddo
3109       return
3110       end subroutine eigtransf1
3111 !-----------------------------------------------------------------------------
3112       subroutine eigtransf(n,ndim,a,b,d,c)
3113
3114 !el      implicit none
3115       integer :: n,ndim
3116       real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3117       integer :: i,j,k
3118       do i=1,n
3119         do j=1,n
3120           c(i,j)=0.0d0
3121           do k=1,n
3122             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3123           enddo
3124         enddo
3125       enddo
3126       return
3127       end subroutine eigtransf
3128 !-----------------------------------------------------------------------------
3129       subroutine sd_verlet1
3130
3131 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3132       use energy_data
3133 !      implicit real*8 (a-h,o-z)
3134 !      include 'DIMENSIONS'
3135 !      include 'COMMON.CONTROL'
3136 !      include 'COMMON.VAR'
3137 !      include 'COMMON.MD'
3138 !#ifndef LANG0
3139 !      include 'COMMON.LANGEVIN'
3140 !#else
3141 !      include 'COMMON.LANGEVIN.lang0'
3142 !#endif
3143 !      include 'COMMON.CHAIN'
3144 !      include 'COMMON.DERIV'
3145 !      include 'COMMON.GEO'
3146 !      include 'COMMON.LOCAL'
3147 !      include 'COMMON.INTERACT'
3148 !      include 'COMMON.IOUNITS'
3149 !      include 'COMMON.NAMES'
3150 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3151 !el      common /stochcalc/ stochforcvec
3152       logical :: lprn = .false.
3153       real(kind=8) :: ddt1,ddt2
3154       integer :: i,j,ind,inres
3155
3156 !      write (iout,*) "dc_old"
3157 !      do i=0,nres
3158 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3159 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3160 !      enddo
3161       do j=1,3
3162         dc_work(j)=dc_old(j,0)
3163         d_t_work(j)=d_t_old(j,0)
3164         d_a_work(j)=d_a_old(j,0)
3165       enddo
3166       ind=3
3167       do i=nnt,nct-1
3168         do j=1,3
3169           dc_work(ind+j)=dc_old(j,i)
3170           d_t_work(ind+j)=d_t_old(j,i)
3171           d_a_work(ind+j)=d_a_old(j,i)
3172         enddo
3173         ind=ind+3
3174       enddo
3175       do i=nnt,nct
3176         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3177           do j=1,3
3178             dc_work(ind+j)=dc_old(j,i+nres)
3179             d_t_work(ind+j)=d_t_old(j,i+nres)
3180             d_a_work(ind+j)=d_a_old(j,i+nres)
3181           enddo
3182           ind=ind+3
3183         endif
3184       enddo
3185 #ifndef LANG0
3186       if (lprn) then
3187       write (iout,*) &
3188         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3189         " vrand_mat2"
3190       do i=1,dimen
3191         do j=1,dimen
3192           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3193             vfric_mat(i,j),afric_mat(i,j),&
3194             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3195         enddo
3196       enddo
3197       endif
3198       do i=1,dimen
3199         ddt1=0.0d0
3200         ddt2=0.0d0
3201         do j=1,dimen
3202           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3203             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3204           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3205           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3206         enddo
3207         d_t_work_new(i)=ddt1+0.5d0*ddt2
3208         d_t_work(i)=ddt1+ddt2
3209       enddo
3210 #endif
3211       do j=1,3
3212         dc(j,0)=dc_work(j)
3213         d_t(j,0)=d_t_work(j)
3214       enddo
3215       ind=3
3216       do i=nnt,nct-1
3217         do j=1,3
3218           dc(j,i)=dc_work(ind+j)
3219           d_t(j,i)=d_t_work(ind+j)
3220         enddo
3221         ind=ind+3
3222       enddo
3223       do i=nnt,nct
3224         if (itype(i).ne.10) then
3225           inres=i+nres
3226           do j=1,3
3227             dc(j,inres)=dc_work(ind+j)
3228             d_t(j,inres)=d_t_work(ind+j)
3229           enddo
3230           ind=ind+3
3231         endif
3232       enddo
3233       return
3234       end subroutine sd_verlet1
3235 !-----------------------------------------------------------------------------
3236       subroutine sd_verlet2
3237
3238 !  Calculating the adjusted velocities for accelerations
3239       use energy_data
3240 !      implicit real*8 (a-h,o-z)
3241 !      include 'DIMENSIONS'
3242 !      include 'COMMON.CONTROL'
3243 !      include 'COMMON.VAR'
3244 !      include 'COMMON.MD'
3245 !#ifndef LANG0
3246 !      include 'COMMON.LANGEVIN'
3247 !#else
3248 !      include 'COMMON.LANGEVIN.lang0'
3249 !#endif
3250 !      include 'COMMON.CHAIN'
3251 !      include 'COMMON.DERIV'
3252 !      include 'COMMON.GEO'
3253 !      include 'COMMON.LOCAL'
3254 !      include 'COMMON.INTERACT'
3255 !      include 'COMMON.IOUNITS'
3256 !      include 'COMMON.NAMES'
3257 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3258        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3259 !el      common /stochcalc/ stochforcvec
3260 !
3261       real(kind=8) :: ddt1,ddt2
3262       integer :: i,j,ind,inres
3263 ! Compute the stochastic forces which contribute to velocity change
3264 !
3265       call stochastic_force(stochforcvecV)
3266
3267 #ifndef LANG0
3268       do i=1,dimen
3269         ddt1=0.0d0
3270         ddt2=0.0d0
3271         do j=1,dimen
3272           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3273           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3274            vrand_mat2(i,j)*stochforcvecV(j)
3275         enddo
3276         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3277       enddo
3278 #endif
3279       do j=1,3
3280         d_t(j,0)=d_t_work(j)
3281       enddo
3282       ind=3
3283       do i=nnt,nct-1
3284         do j=1,3
3285           d_t(j,i)=d_t_work(ind+j)
3286         enddo
3287         ind=ind+3
3288       enddo
3289       do i=nnt,nct
3290         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3291           inres=i+nres
3292           do j=1,3
3293             d_t(j,inres)=d_t_work(ind+j)
3294           enddo
3295           ind=ind+3
3296         endif
3297       enddo
3298       return
3299       end subroutine sd_verlet2
3300 !-----------------------------------------------------------------------------
3301       subroutine sd_verlet_ciccotti_setup
3302
3303 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3304 ! version
3305 !      implicit real*8 (a-h,o-z)
3306 !      include 'DIMENSIONS'
3307       use control, only: tcpu
3308       use control_data
3309 #ifdef MPI
3310       include 'mpif.h'
3311 #endif
3312 !      include 'COMMON.CONTROL'
3313 !      include 'COMMON.VAR'
3314 !      include 'COMMON.MD'
3315 !#ifndef LANG0
3316 !      include 'COMMON.LANGEVIN'
3317 !#else
3318 !      include 'COMMON.LANGEVIN.lang0'
3319 !#endif
3320 !      include 'COMMON.CHAIN'
3321 !      include 'COMMON.DERIV'
3322 !      include 'COMMON.GEO'
3323 !      include 'COMMON.LOCAL'
3324 !      include 'COMMON.INTERACT'
3325 !      include 'COMMON.IOUNITS'
3326 !      include 'COMMON.NAMES'
3327 !      include 'COMMON.TIME1'
3328       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
3329       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3330       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3331         prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3332       logical :: lprn = .false.
3333       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3334       real(kind=8) :: ktm,gdt,egdt,tt0
3335       integer :: i,maxres2
3336 #ifdef MPI
3337       tt0 = MPI_Wtime()
3338 #else
3339       tt0 = tcpu()
3340 #endif
3341 !
3342 ! AL 8/17/04 Code adapted from tinker
3343 !
3344 ! Get the frictional and random terms for stochastic dynamics in the
3345 ! eigenspace of mass-scaled UNRES friction matrix
3346 !
3347       maxres2=2*nres
3348       do i = 1, dimen
3349             write (iout,*) "i",i," fricgam",fricgam(i)
3350             gdt = fricgam(i) * d_time
3351 !
3352 ! Stochastic dynamics reduces to simple MD for zero friction
3353 !
3354             if (gdt .le. zero) then
3355                pfric_vec(i) = 1.0d0
3356                vfric_vec(i) = d_time
3357                afric_vec(i) = 0.5d0*d_time*d_time
3358                prand_vec(i) = afric_vec(i)
3359                vrand_vec2(i) = vfric_vec(i)
3360 !
3361 ! Analytical expressions when friction coefficient is large
3362 !
3363             else
3364                egdt = dexp(-gdt)
3365                pfric_vec(i) = egdt
3366                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3367                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3368                prand_vec(i) = afric_vec(i)
3369                vrand_vec2(i) = vfric_vec(i)
3370 !
3371 ! Compute the scaling factors of random terms for the nonzero friction case
3372 !
3373 !               ktm = 0.5d0*d_time/fricgam(i)
3374 !               psig = dsqrt(ktm*pterm) / fricgam(i)
3375 !               vsig = dsqrt(ktm*vterm)
3376 !               prand_vec(i) = psig*afric_vec(i)
3377 !               vrand_vec2(i) = vsig*vfric_vec(i)
3378             end if
3379       end do
3380       if (lprn) then
3381       write (iout,*) &
3382         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3383         " vrand_vec2"
3384       do i=1,dimen
3385         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3386             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3387       enddo
3388       endif
3389 !
3390 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3391 !
3392       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3393       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3394       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3395       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3396       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3397 #ifdef MPI
3398       t_sdsetup=t_sdsetup+MPI_Wtime()
3399 #else
3400       t_sdsetup=t_sdsetup+tcpu()-tt0
3401 #endif
3402       return
3403       end subroutine sd_verlet_ciccotti_setup
3404 !-----------------------------------------------------------------------------
3405       subroutine sd_verlet1_ciccotti
3406
3407 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3408 !      implicit real*8 (a-h,o-z)
3409       use energy_data
3410 !      include 'DIMENSIONS'
3411 #ifdef MPI
3412       include 'mpif.h'
3413 #endif
3414 !      include 'COMMON.CONTROL'
3415 !      include 'COMMON.VAR'
3416 !      include 'COMMON.MD'
3417 !#ifndef LANG0
3418 !      include 'COMMON.LANGEVIN'
3419 !#else
3420 !      include 'COMMON.LANGEVIN.lang0'
3421 !#endif
3422 !      include 'COMMON.CHAIN'
3423 !      include 'COMMON.DERIV'
3424 !      include 'COMMON.GEO'
3425 !      include 'COMMON.LOCAL'
3426 !      include 'COMMON.INTERACT'
3427 !      include 'COMMON.IOUNITS'
3428 !      include 'COMMON.NAMES'
3429 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3430 !el      common /stochcalc/ stochforcvec
3431       logical :: lprn = .false.
3432       real(kind=8) :: ddt1,ddt2
3433       integer :: i,j,ind,inres
3434 !      write (iout,*) "dc_old"
3435 !      do i=0,nres
3436 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3437 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3438 !      enddo
3439       do j=1,3
3440         dc_work(j)=dc_old(j,0)
3441         d_t_work(j)=d_t_old(j,0)
3442         d_a_work(j)=d_a_old(j,0)
3443       enddo
3444       ind=3
3445       do i=nnt,nct-1
3446         do j=1,3
3447           dc_work(ind+j)=dc_old(j,i)
3448           d_t_work(ind+j)=d_t_old(j,i)
3449           d_a_work(ind+j)=d_a_old(j,i)
3450         enddo
3451         ind=ind+3
3452       enddo
3453       do i=nnt,nct
3454         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3455           do j=1,3
3456             dc_work(ind+j)=dc_old(j,i+nres)
3457             d_t_work(ind+j)=d_t_old(j,i+nres)
3458             d_a_work(ind+j)=d_a_old(j,i+nres)
3459           enddo
3460           ind=ind+3
3461         endif
3462       enddo
3463
3464 #ifndef LANG0
3465       if (lprn) then
3466       write (iout,*) &
3467         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3468         " vrand_mat2"
3469       do i=1,dimen
3470         do j=1,dimen
3471           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3472             vfric_mat(i,j),afric_mat(i,j),&
3473             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3474         enddo
3475       enddo
3476       endif
3477       do i=1,dimen
3478         ddt1=0.0d0
3479         ddt2=0.0d0
3480         do j=1,dimen
3481           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3482             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3483           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3484           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3485         enddo
3486         d_t_work_new(i)=ddt1+0.5d0*ddt2
3487         d_t_work(i)=ddt1+ddt2
3488       enddo
3489 #endif
3490       do j=1,3
3491         dc(j,0)=dc_work(j)
3492         d_t(j,0)=d_t_work(j)
3493       enddo
3494       ind=3
3495       do i=nnt,nct-1
3496         do j=1,3
3497           dc(j,i)=dc_work(ind+j)
3498           d_t(j,i)=d_t_work(ind+j)
3499         enddo
3500         ind=ind+3
3501       enddo
3502       do i=nnt,nct
3503         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3504           inres=i+nres
3505           do j=1,3
3506             dc(j,inres)=dc_work(ind+j)
3507             d_t(j,inres)=d_t_work(ind+j)
3508           enddo
3509           ind=ind+3
3510         endif
3511       enddo
3512       return
3513       end subroutine sd_verlet1_ciccotti
3514 !-----------------------------------------------------------------------------
3515       subroutine sd_verlet2_ciccotti
3516
3517 !  Calculating the adjusted velocities for accelerations
3518       use energy_data
3519 !      implicit real*8 (a-h,o-z)
3520 !      include 'DIMENSIONS'
3521 !      include 'COMMON.CONTROL'
3522 !      include 'COMMON.VAR'
3523 !      include 'COMMON.MD'
3524 !#ifndef LANG0
3525 !      include 'COMMON.LANGEVIN'
3526 !#else
3527 !      include 'COMMON.LANGEVIN.lang0'
3528 !#endif
3529 !      include 'COMMON.CHAIN'
3530 !      include 'COMMON.DERIV'
3531 !      include 'COMMON.GEO'
3532 !      include 'COMMON.LOCAL'
3533 !      include 'COMMON.INTERACT'
3534 !      include 'COMMON.IOUNITS'
3535 !      include 'COMMON.NAMES'
3536 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3537        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3538 !el      common /stochcalc/ stochforcvec
3539       real(kind=8) :: ddt1,ddt2
3540       integer :: i,j,ind,inres
3541 !
3542 ! Compute the stochastic forces which contribute to velocity change
3543 !
3544       call stochastic_force(stochforcvecV)
3545 #ifndef LANG0
3546       do i=1,dimen
3547         ddt1=0.0d0
3548         ddt2=0.0d0
3549         do j=1,dimen
3550
3551           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3552 !          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3553           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3554         enddo
3555         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3556       enddo
3557 #endif
3558       do j=1,3
3559         d_t(j,0)=d_t_work(j)
3560       enddo
3561       ind=3
3562       do i=nnt,nct-1
3563         do j=1,3
3564           d_t(j,i)=d_t_work(ind+j)
3565         enddo
3566         ind=ind+3
3567       enddo
3568       do i=nnt,nct
3569         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3570           inres=i+nres
3571           do j=1,3
3572             d_t(j,inres)=d_t_work(ind+j)
3573           enddo
3574           ind=ind+3
3575         endif
3576       enddo
3577       return
3578       end subroutine sd_verlet2_ciccotti
3579 #endif
3580 !-----------------------------------------------------------------------------
3581 ! moments.f
3582 !-----------------------------------------------------------------------------
3583       subroutine inertia_tensor
3584
3585 ! Calculating the intertia tensor for the entire protein in order to
3586 ! remove the perpendicular components of velocity matrix which cause
3587 ! the molecule to rotate.
3588       use comm_gucio
3589       use energy_data
3590 !       implicit real*8 (a-h,o-z)
3591 !       include 'DIMENSIONS'
3592 !       include 'COMMON.CONTROL'
3593 !       include 'COMMON.VAR'
3594 !       include 'COMMON.MD'
3595 !       include 'COMMON.CHAIN'
3596 !       include 'COMMON.DERIV'
3597 !       include 'COMMON.GEO'
3598 !       include 'COMMON.LOCAL'
3599 !       include 'COMMON.INTERACT'
3600 !       include 'COMMON.IOUNITS'
3601 !       include 'COMMON.NAMES'
3602
3603       real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3604       real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3605       real(kind=8) :: M_SC,mag,mag2
3606       real(kind=8),dimension(3,0:nres) :: vpp   !(3,0:MAXRES)
3607       real(kind=8),dimension(3) :: vs_p,pp,incr,v
3608       real(kind=8),dimension(3,3) :: pr1,pr2
3609
3610 !el      common /gucio/ cm
3611       integer :: iti,inres,i,j,k
3612         do i=1,3
3613            do j=1,3
3614              Im(i,j)=0.0d0
3615              pr1(i,j)=0.0d0
3616              pr2(i,j)=0.0d0
3617            enddo
3618            L(i)=0.0d0
3619            cm(i)=0.0d0
3620            vrot(i)=0.0d0
3621         enddo
3622 !   calculating the center of the mass of the protein
3623         do i=nnt,nct-1
3624           do j=1,3
3625             cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
3626           enddo
3627         enddo
3628         do j=1,3
3629          cm(j)=mp*cm(j)
3630         enddo
3631         M_SC=0.0d0
3632         do i=nnt,nct
3633            iti=iabs(itype(i))
3634            M_SC=M_SC+msc(iabs(iti))
3635            inres=i+nres
3636            do j=1,3
3637             cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
3638            enddo
3639         enddo
3640         do j=1,3
3641           cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
3642         enddo
3643
3644         do i=nnt,nct-1
3645           do j=1,3
3646             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3647           enddo
3648           Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
3649           Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
3650           Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
3651           Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
3652           Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
3653           Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
3654         enddo
3655
3656         do i=nnt,nct
3657            iti=iabs(itype(i))
3658            inres=i+nres
3659            do j=1,3
3660              pr(j)=c(j,inres)-cm(j)
3661            enddo
3662           Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
3663           Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
3664           Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
3665           Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
3666           Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
3667           Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
3668         enddo
3669
3670         do i=nnt,nct-1
3671           Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
3672           vbld(i+1)*vbld(i+1)*0.25d0
3673           Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* &
3674           vbld(i+1)*vbld(i+1)*0.25d0
3675           Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* &
3676           vbld(i+1)*vbld(i+1)*0.25d0
3677           Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* &
3678           vbld(i+1)*vbld(i+1)*0.25d0
3679           Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* &
3680           vbld(i+1)*vbld(i+1)*0.25d0
3681           Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* &
3682           vbld(i+1)*vbld(i+1)*0.25d0
3683         enddo
3684
3685
3686         do i=nnt,nct
3687          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3688            iti=iabs(itype(i))
3689            inres=i+nres
3690           Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
3691           dc_norm(1,inres))*vbld(inres)*vbld(inres)
3692           Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* &
3693           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3694           Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
3695           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3696           Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
3697           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3698           Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* &
3699           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3700           Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* &
3701           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3702          endif
3703         enddo
3704
3705         call angmom(cm,L)
3706 !        write(iout,*) "The angular momentum before adjustment:"
3707 !        write(iout,*) (L(j),j=1,3)
3708
3709         Im(2,1)=Im(1,2)
3710         Im(3,1)=Im(1,3)
3711         Im(3,2)=Im(2,3)
3712
3713 !  Copying the Im matrix for the djacob subroutine
3714         do i=1,3
3715           do j=1,3
3716             Imcp(i,j)=Im(i,j)
3717             Id(i,j)=0.0d0
3718           enddo
3719         enddo
3720
3721 !   Finding the eigenvectors and eignvalues of the inertia tensor
3722        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3723 !       write (iout,*) "Eigenvalues & Eigenvectors"
3724 !       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3725 !       write (iout,*)
3726 !       do i=1,3
3727 !         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3728 !       enddo
3729 !   Constructing the diagonalized matrix
3730        do i=1,3
3731          if (dabs(eigval(i)).gt.1.0d-15) then
3732            Id(i,i)=1.0d0/eigval(i)
3733          else
3734            Id(i,i)=0.0d0
3735          endif
3736        enddo
3737         do i=1,3
3738            do j=1,3
3739               Imcp(i,j)=eigvec(j,i)
3740            enddo
3741         enddo
3742         do i=1,3
3743            do j=1,3
3744               do k=1,3
3745                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3746               enddo
3747            enddo
3748         enddo
3749         do i=1,3
3750            do j=1,3
3751               do k=1,3
3752                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3753               enddo
3754            enddo
3755         enddo
3756 !  Calculating the total rotational velocity of the molecule
3757        do i=1,3
3758          do j=1,3
3759            vrot(i)=vrot(i)+pr2(i,j)*L(j)
3760          enddo
3761        enddo
3762 !   Resetting the velocities
3763        do i=nnt,nct-1
3764          call vecpr(vrot(1),dc(1,i),vp)
3765          do j=1,3
3766            d_t(j,i)=d_t(j,i)-vp(j)
3767           enddo
3768         enddo
3769         do i=nnt,nct
3770         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3771            inres=i+nres
3772            call vecpr(vrot(1),dc(1,inres),vp)
3773            do j=1,3
3774              d_t(j,inres)=d_t(j,inres)-vp(j)
3775            enddo
3776         endif
3777        enddo
3778        call angmom(cm,L)
3779 !       write(iout,*) "The angular momentum after adjustment:"
3780 !       write(iout,*) (L(j),j=1,3)
3781
3782       return
3783       end subroutine inertia_tensor
3784 !-----------------------------------------------------------------------------
3785       subroutine angmom(cm,L)
3786
3787       use energy_data
3788 !       implicit real*8 (a-h,o-z)
3789 !       include 'DIMENSIONS'
3790 !       include 'COMMON.CONTROL'
3791 !       include 'COMMON.VAR'
3792 !       include 'COMMON.MD'
3793 !       include 'COMMON.CHAIN'
3794 !       include 'COMMON.DERIV'
3795 !       include 'COMMON.GEO'
3796 !       include 'COMMON.LOCAL'
3797 !       include 'COMMON.INTERACT'
3798 !       include 'COMMON.IOUNITS'
3799 !       include 'COMMON.NAMES'
3800
3801       real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3802       integer :: iti,inres,i,j
3803 !  Calculate the angular momentum
3804        do j=1,3
3805           L(j)=0.0d0
3806        enddo
3807        do j=1,3
3808           incr(j)=d_t(j,0)
3809        enddo
3810        do i=nnt,nct-1
3811           do j=1,3
3812             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3813           enddo
3814           do j=1,3
3815             v(j)=incr(j)+0.5d0*d_t(j,i)
3816           enddo
3817           do j=1,3
3818             incr(j)=incr(j)+d_t(j,i)
3819           enddo
3820           call vecpr(pr(1),v(1),vp)
3821           do j=1,3
3822             L(j)=L(j)+mp*vp(j)
3823           enddo
3824           do j=1,3
3825              pr(j)=0.5d0*dc(j,i)
3826              pp(j)=0.5d0*d_t(j,i)
3827           enddo
3828          call vecpr(pr(1),pp(1),vp)
3829          do j=1,3
3830              L(j)=L(j)+Ip*vp(j)
3831           enddo
3832         enddo
3833         do j=1,3
3834           incr(j)=d_t(j,0)
3835         enddo
3836         do i=nnt,nct
3837          iti=iabs(itype(i))
3838          inres=i+nres
3839          do j=1,3
3840            pr(j)=c(j,inres)-cm(j)
3841          enddo
3842          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3843            do j=1,3
3844              v(j)=incr(j)+d_t(j,inres)
3845            enddo
3846          else
3847            do j=1,3
3848              v(j)=incr(j)
3849            enddo
3850          endif
3851          call vecpr(pr(1),v(1),vp)
3852 !         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3853 !     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3854          do j=1,3
3855             L(j)=L(j)+msc(iabs(iti))*vp(j)
3856          enddo
3857 !         write (iout,*) "L",(l(j),j=1,3)
3858          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3859            do j=1,3
3860             v(j)=incr(j)+d_t(j,inres)
3861            enddo
3862            call vecpr(dc(1,inres),d_t(1,inres),vp)
3863            do j=1,3
3864              L(j)=L(j)+Isc(iti)*vp(j)
3865           enddo
3866          endif
3867          do j=1,3
3868              incr(j)=incr(j)+d_t(j,i)
3869          enddo
3870        enddo
3871       return
3872       end subroutine angmom
3873 !-----------------------------------------------------------------------------
3874       subroutine vcm_vel(vcm)
3875
3876       use energy_data
3877 !       implicit real*8 (a-h,o-z)
3878 !       include 'DIMENSIONS'
3879 !       include 'COMMON.VAR'
3880 !       include 'COMMON.MD'
3881 !       include 'COMMON.CHAIN'
3882 !       include 'COMMON.DERIV'
3883 !       include 'COMMON.GEO'
3884 !       include 'COMMON.LOCAL'
3885 !       include 'COMMON.INTERACT'
3886 !       include 'COMMON.IOUNITS'
3887        real(kind=8),dimension(3) :: vcm,vv
3888        real(kind=8) :: summas,amas
3889        integer :: i,j
3890
3891        do j=1,3
3892          vcm(j)=0.0d0
3893          vv(j)=d_t(j,0)
3894        enddo
3895        summas=0.0d0
3896        do i=nnt,nct
3897          if (i.lt.nct) then
3898            summas=summas+mp
3899            do j=1,3
3900              vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
3901            enddo
3902          endif
3903          amas=msc(iabs(itype(i)))
3904          summas=summas+amas
3905          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3906            do j=1,3
3907              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
3908            enddo
3909          else
3910            do j=1,3
3911              vcm(j)=vcm(j)+amas*vv(j)
3912            enddo
3913          endif
3914          do j=1,3
3915            vv(j)=vv(j)+d_t(j,i)
3916          enddo
3917        enddo
3918 !       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
3919        do j=1,3
3920          vcm(j)=vcm(j)/summas
3921        enddo
3922       return
3923       end subroutine vcm_vel
3924 !-----------------------------------------------------------------------------
3925 ! rattle.F
3926 !-----------------------------------------------------------------------------
3927       subroutine rattle1
3928 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
3929 ! AL 9/24/04
3930       use comm_przech
3931       use energy_data
3932 !      implicit real*8 (a-h,o-z)
3933 !      include 'DIMENSIONS'
3934 #ifdef RATTLE
3935 !      include 'COMMON.CONTROL'
3936 !      include 'COMMON.VAR'
3937 !      include 'COMMON.MD'
3938 !#ifndef LANG0
3939 !      include 'COMMON.LANGEVIN'
3940 !#else
3941 !      include 'COMMON.LANGEVIN.lang0'
3942 !#endif
3943 !      include 'COMMON.CHAIN'
3944 !      include 'COMMON.DERIV'
3945 !      include 'COMMON.GEO'
3946 !      include 'COMMON.LOCAL'
3947 !      include 'COMMON.INTERACT'
3948 !      include 'COMMON.IOUNITS'
3949 !      include 'COMMON.NAMES'
3950 !      include 'COMMON.TIME1'
3951 !el      real(kind=8) :: gginv(2*nres,2*nres),&
3952 !el       gdc(3,2*nres,2*nres)
3953       real(kind=8) :: dC_uncor(3,2*nres)        !,&
3954 !el      real(kind=8) :: Cmat(2*nres,2*nres)
3955       real(kind=8) :: x(2*nres),xcorr(3,2*nres)         !maxres2=2*maxres
3956 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
3957 !el      common /przechowalnia/ nbond
3958       integer :: max_rattle = 5
3959       logical :: lprn = .false., lprn1 = .false., not_done
3960       real(kind=8) :: tol_rattle = 1.0d-5
3961
3962       integer :: ii,i,j,jj,l,ind,ind1,nres2
3963       nres2=2*nres
3964
3965 !el /common/ przechowalnia
3966
3967       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
3968       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
3969       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
3970 !el--------
3971       if (lprn) write (iout,*) "RATTLE1"
3972       nbond=nct-nnt
3973       do i=nnt,nct
3974         if (itype(i).ne.10) nbond=nbond+1
3975       enddo
3976 ! Make a folded form of the Ginv-matrix
3977       ind=0
3978       ii=0
3979       do i=nnt,nct-1
3980         ii=ii+1
3981         do j=1,3
3982           ind=ind+1
3983           ind1=0
3984           jj=0
3985           do k=nnt,nct-1
3986             jj=jj+1
3987             do l=1,3
3988               ind1=ind1+1
3989               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3990             enddo
3991           enddo
3992           do k=nnt,nct
3993             if (itype(k).ne.10) then
3994               jj=jj+1
3995               do l=1,3
3996                 ind1=ind1+1
3997                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3998               enddo
3999             endif
4000           enddo
4001         enddo
4002       enddo
4003       do i=nnt,nct
4004         if (itype(i).ne.10) then
4005           ii=ii+1
4006           do j=1,3
4007             ind=ind+1
4008             ind1=0
4009             jj=0
4010             do k=nnt,nct-1
4011               jj=jj+1
4012               do l=1,3
4013                 ind1=ind1+1
4014                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4015               enddo
4016             enddo
4017             do k=nnt,nct
4018               if (itype(k).ne.10) then
4019                 jj=jj+1
4020                 do l=1,3
4021                   ind1=ind1+1
4022                   if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4023                 enddo
4024               endif
4025             enddo
4026           enddo
4027         endif
4028       enddo
4029       if (lprn1) then
4030         write (iout,*) "Matrix GGinv"
4031         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4032       endif
4033       not_done=.true.
4034       iter=0
4035       do while (not_done)
4036       iter=iter+1
4037       if (iter.gt.max_rattle) then
4038         write (iout,*) "Error - too many iterations in RATTLE."
4039         stop
4040       endif
4041 ! Calculate the matrix C = GG**(-1) dC_old o dC
4042       ind1=0
4043       do i=nnt,nct-1
4044         ind1=ind1+1
4045         do j=1,3
4046           dC_uncor(j,ind1)=dC(j,i)
4047         enddo
4048       enddo
4049       do i=nnt,nct
4050         if (itype(i).ne.10) then
4051           ind1=ind1+1
4052           do j=1,3
4053             dC_uncor(j,ind1)=dC(j,i+nres)
4054           enddo
4055         endif
4056       enddo
4057       do i=1,nbond
4058         ind=0
4059         do k=nnt,nct-1
4060           ind=ind+1
4061           do j=1,3
4062             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4063           enddo
4064         enddo
4065         do k=nnt,nct
4066           if (itype(k).ne.10) then
4067             ind=ind+1
4068             do j=1,3
4069               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4070             enddo
4071           endif
4072         enddo
4073       enddo
4074 ! Calculate deviations from standard virtual-bond lengths
4075       ind=0
4076       do i=nnt,nct-1
4077         ind=ind+1
4078         x(ind)=vbld(i+1)**2-vbl**2
4079       enddo
4080       do i=nnt,nct
4081         if (itype(i).ne.10) then
4082           ind=ind+1
4083           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4084         endif
4085       enddo
4086       if (lprn) then
4087         write (iout,*) "Coordinates and violations"
4088         do i=1,nbond
4089           write(iout,'(i5,3f10.5,5x,e15.5)') &
4090            i,(dC_uncor(j,i),j=1,3),x(i)
4091         enddo
4092         write (iout,*) "Velocities and violations"
4093         ind=0
4094         do i=nnt,nct-1
4095           ind=ind+1
4096           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4097            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4098         enddo
4099         do i=nnt,nct
4100           if (itype(i).ne.10) then
4101             ind=ind+1
4102             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4103              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4104              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4105           endif
4106         enddo
4107 !        write (iout,*) "gdc"
4108 !        do i=1,nbond
4109 !          write (iout,*) "i",i
4110 !          do j=1,nbond
4111 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4112 !          enddo
4113 !        enddo
4114       endif
4115       xmax=dabs(x(1))
4116       do i=2,nbond
4117         if (dabs(x(i)).gt.xmax) then
4118           xmax=dabs(x(i))
4119         endif
4120       enddo
4121       if (xmax.lt.tol_rattle) then
4122         not_done=.false.
4123         goto 100
4124       endif
4125 ! Calculate the matrix of the system of equations
4126       do i=1,nbond
4127         do j=1,nbond
4128           Cmat(i,j)=0.0d0
4129           do k=1,3
4130             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4131           enddo
4132         enddo
4133       enddo
4134       if (lprn1) then
4135         write (iout,*) "Matrix Cmat"
4136         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4137       endif
4138       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4139 ! Add constraint term to positions
4140       ind=0
4141       do i=nnt,nct-1
4142         ind=ind+1
4143         do j=1,3
4144           xx=0.0d0
4145           do ii=1,nbond
4146             xx = xx+x(ii)*gdc(j,ind,ii)
4147           enddo
4148           xx=0.5d0*xx
4149           dC(j,i)=dC(j,i)-xx
4150           d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4151         enddo
4152       enddo
4153       do i=nnt,nct
4154         if (itype(i).ne.10) then
4155           ind=ind+1
4156           do j=1,3
4157             xx=0.0d0
4158             do ii=1,nbond
4159               xx = xx+x(ii)*gdc(j,ind,ii)
4160             enddo
4161             xx=0.5d0*xx
4162             dC(j,i+nres)=dC(j,i+nres)-xx
4163             d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4164           enddo
4165         endif
4166       enddo
4167 ! Rebuild the chain using the new coordinates
4168       call chainbuild_cart
4169       if (lprn) then
4170         write (iout,*) "New coordinates, Lagrange multipliers,",&
4171         " and differences between actual and standard bond lengths"
4172         ind=0
4173         do i=nnt,nct-1
4174           ind=ind+1
4175           xx=vbld(i+1)**2-vbl**2
4176           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4177               i,(dC(j,i),j=1,3),x(ind),xx
4178         enddo
4179         do i=nnt,nct
4180           if (itype(i).ne.10) then
4181             ind=ind+1
4182             xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4183             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4184              i,(dC(j,i+nres),j=1,3),x(ind),xx
4185           endif
4186         enddo
4187         write (iout,*) "Velocities and violations"
4188         ind=0
4189         do i=nnt,nct-1
4190           ind=ind+1
4191           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4192            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4193         enddo
4194         do i=nnt,nct
4195           if (itype(i).ne.10) then
4196             ind=ind+1
4197             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4198              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4199              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4200           endif
4201         enddo
4202       endif
4203       enddo
4204   100 continue
4205       return
4206    10 write (iout,*) "Error - singularity in solving the system",&
4207        " of equations for Lagrange multipliers."
4208       stop
4209 #else
4210       write (iout,*) &
4211        "RATTLE inactive; use -DRATTLE switch at compile time."
4212       stop
4213 #endif
4214       end subroutine rattle1
4215 !-----------------------------------------------------------------------------
4216       subroutine rattle2
4217 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4218 ! AL 9/24/04
4219       use comm_przech
4220       use energy_data
4221 !      implicit real*8 (a-h,o-z)
4222 !      include 'DIMENSIONS'
4223 #ifdef RATTLE
4224 !      include 'COMMON.CONTROL'
4225 !      include 'COMMON.VAR'
4226 !      include 'COMMON.MD'
4227 !#ifndef LANG0
4228 !      include 'COMMON.LANGEVIN'
4229 !#else
4230 !      include 'COMMON.LANGEVIN.lang0'
4231 !#endif
4232 !      include 'COMMON.CHAIN'
4233 !      include 'COMMON.DERIV'
4234 !      include 'COMMON.GEO'
4235 !      include 'COMMON.LOCAL'
4236 !      include 'COMMON.INTERACT'
4237 !      include 'COMMON.IOUNITS'
4238 !      include 'COMMON.NAMES'
4239 !      include 'COMMON.TIME1'
4240 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4241 !el       gdc(3,2*nres,2*nres)
4242       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4243 !el       Cmat(2*nres,2*nres)
4244       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4245 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4246 !el      common /przechowalnia/ nbond
4247       integer :: max_rattle = 5
4248       logical :: lprn = .false., lprn1 = .false., not_done
4249       real(kind=8) :: tol_rattle = 1.0d-5
4250       integer :: nres2
4251       nres2=2*nres
4252
4253 !el /common/ przechowalnia
4254       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4255       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4256       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4257 !el--------
4258       if (lprn) write (iout,*) "RATTLE2"
4259       if (lprn) write (iout,*) "Velocity correction"
4260 ! Calculate the matrix G dC
4261       do i=1,nbond
4262         ind=0
4263         do k=nnt,nct-1
4264           ind=ind+1
4265           do j=1,3
4266             gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4267           enddo
4268         enddo
4269         do k=nnt,nct
4270           if (itype(k).ne.10) then
4271             ind=ind+1
4272             do j=1,3
4273               gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4274             enddo
4275           endif
4276         enddo
4277       enddo
4278 !      if (lprn) then
4279 !        write (iout,*) "gdc"
4280 !        do i=1,nbond
4281 !          write (iout,*) "i",i
4282 !          do j=1,nbond
4283 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4284 !          enddo
4285 !        enddo
4286 !      endif
4287 ! Calculate the matrix of the system of equations
4288       ind=0
4289       do i=nnt,nct-1
4290         ind=ind+1
4291         do j=1,nbond
4292           Cmat(ind,j)=0.0d0
4293           do k=1,3
4294             Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4295           enddo
4296         enddo
4297       enddo
4298       do i=nnt,nct
4299         if (itype(i).ne.10) then
4300           ind=ind+1
4301           do j=1,nbond
4302             Cmat(ind,j)=0.0d0
4303             do k=1,3
4304               Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4305             enddo
4306           enddo
4307         endif
4308       enddo
4309 ! Calculate the scalar product dC o d_t_new
4310       ind=0
4311       do i=nnt,nct-1
4312         ind=ind+1
4313         x(ind)=scalar(d_t(1,i),dC(1,i))
4314       enddo
4315       do i=nnt,nct
4316         if (itype(i).ne.10) then
4317           ind=ind+1
4318           x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4319         endif
4320       enddo
4321       if (lprn) then
4322         write (iout,*) "Velocities and violations"
4323         ind=0
4324         do i=nnt,nct-1
4325           ind=ind+1
4326           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4327            i,ind,(d_t(j,i),j=1,3),x(ind)
4328         enddo
4329         do i=nnt,nct
4330           if (itype(i).ne.10) then
4331             ind=ind+1
4332             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4333              i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4334           endif
4335         enddo
4336       endif
4337       xmax=dabs(x(1))
4338       do i=2,nbond
4339         if (dabs(x(i)).gt.xmax) then
4340           xmax=dabs(x(i))
4341         endif
4342       enddo
4343       if (xmax.lt.tol_rattle) then
4344         not_done=.false.
4345         goto 100
4346       endif
4347       if (lprn1) then
4348         write (iout,*) "Matrix Cmat"
4349         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4350       endif
4351       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4352 ! Add constraint term to velocities
4353       ind=0
4354       do i=nnt,nct-1
4355         ind=ind+1
4356         do j=1,3
4357           xx=0.0d0
4358           do ii=1,nbond
4359             xx = xx+x(ii)*gdc(j,ind,ii)
4360           enddo
4361           d_t(j,i)=d_t(j,i)-xx
4362         enddo
4363       enddo
4364       do i=nnt,nct
4365         if (itype(i).ne.10) then
4366           ind=ind+1
4367           do j=1,3
4368             xx=0.0d0
4369             do ii=1,nbond
4370               xx = xx+x(ii)*gdc(j,ind,ii)
4371             enddo
4372             d_t(j,i+nres)=d_t(j,i+nres)-xx
4373           enddo
4374         endif
4375       enddo
4376       if (lprn) then
4377         write (iout,*) &
4378           "New velocities, Lagrange multipliers violations"
4379         ind=0
4380         do i=nnt,nct-1
4381           ind=ind+1
4382           if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4383              i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4384         enddo
4385         do i=nnt,nct
4386           if (itype(i).ne.10) then
4387             ind=ind+1
4388             write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4389               i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4390               scalar(d_t(1,i+nres),dC(1,i+nres))
4391           endif
4392         enddo
4393       endif
4394   100 continue
4395       return
4396    10 write (iout,*) "Error - singularity in solving the system",&
4397        " of equations for Lagrange multipliers."
4398       stop
4399 #else
4400       write (iout,*) &
4401        "RATTLE inactive; use -DRATTLE option at compile time."
4402       stop
4403 #endif
4404       end subroutine rattle2
4405 !-----------------------------------------------------------------------------
4406       subroutine rattle_brown
4407 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4408 ! AL 9/24/04
4409       use comm_przech
4410       use energy_data
4411 !      implicit real*8 (a-h,o-z)
4412 !      include 'DIMENSIONS'
4413 #ifdef RATTLE
4414 !      include 'COMMON.CONTROL'
4415 !      include 'COMMON.VAR'
4416 !      include 'COMMON.MD'
4417 !#ifndef LANG0
4418 !      include 'COMMON.LANGEVIN'
4419 !#else
4420 !      include 'COMMON.LANGEVIN.lang0'
4421 !#endif
4422 !      include 'COMMON.CHAIN'
4423 !      include 'COMMON.DERIV'
4424 !      include 'COMMON.GEO'
4425 !      include 'COMMON.LOCAL'
4426 !      include 'COMMON.INTERACT'
4427 !      include 'COMMON.IOUNITS'
4428 !      include 'COMMON.NAMES'
4429 !      include 'COMMON.TIME1'
4430 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4431 !el       gdc(3,2*nres,2*nres)
4432       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4433 !el      real(kind=8) :: Cmat(2*nres,2*nres)
4434       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4435 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4436 !el      common /przechowalnia/ nbond
4437       integer :: max_rattle = 5
4438       logical :: lprn = .true., lprn1 = .true., not_done
4439       real(kind=8) :: tol_rattle = 1.0d-5
4440       integer :: nres2
4441       nres2=2*nres
4442
4443 !el /common/ przechowalnia
4444       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4445       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4446       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4447 !el--------
4448
4449       if (lprn) write (iout,*) "RATTLE_BROWN"
4450       nbond=nct-nnt
4451       do i=nnt,nct
4452         if (itype(i).ne.10) nbond=nbond+1
4453       enddo
4454 ! Make a folded form of the Ginv-matrix
4455       ind=0
4456       ii=0
4457       do i=nnt,nct-1
4458         ii=ii+1
4459         do j=1,3
4460           ind=ind+1
4461           ind1=0
4462           jj=0
4463           do k=nnt,nct-1
4464             jj=jj+1
4465             do l=1,3
4466               ind1=ind1+1
4467               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4468             enddo
4469           enddo
4470           do k=nnt,nct
4471             if (itype(k).ne.10) then
4472               jj=jj+1
4473               do l=1,3
4474                 ind1=ind1+1
4475                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4476               enddo
4477             endif
4478           enddo
4479         enddo
4480       enddo
4481       do i=nnt,nct
4482         if (itype(i).ne.10) then
4483           ii=ii+1
4484           do j=1,3
4485             ind=ind+1
4486             ind1=0
4487             jj=0
4488             do k=nnt,nct-1
4489               jj=jj+1
4490               do l=1,3
4491                 ind1=ind1+1
4492                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4493               enddo
4494             enddo
4495             do k=nnt,nct
4496               if (itype(k).ne.10) then
4497                 jj=jj+1
4498                 do l=1,3
4499                   ind1=ind1+1
4500                   if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4501                 enddo
4502               endif
4503             enddo
4504           enddo
4505         endif
4506       enddo
4507       if (lprn1) then
4508         write (iout,*) "Matrix GGinv"
4509         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4510       endif
4511       not_done=.true.
4512       iter=0
4513       do while (not_done)
4514       iter=iter+1
4515       if (iter.gt.max_rattle) then
4516         write (iout,*) "Error - too many iterations in RATTLE."
4517         stop
4518       endif
4519 ! Calculate the matrix C = GG**(-1) dC_old o dC
4520       ind1=0
4521       do i=nnt,nct-1
4522         ind1=ind1+1
4523         do j=1,3
4524           dC_uncor(j,ind1)=dC(j,i)
4525         enddo
4526       enddo
4527       do i=nnt,nct
4528         if (itype(i).ne.10) then
4529           ind1=ind1+1
4530           do j=1,3
4531             dC_uncor(j,ind1)=dC(j,i+nres)
4532           enddo
4533         endif
4534       enddo
4535       do i=1,nbond
4536         ind=0
4537         do k=nnt,nct-1
4538           ind=ind+1
4539           do j=1,3
4540             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4541           enddo
4542         enddo
4543         do k=nnt,nct
4544           if (itype(k).ne.10) then
4545             ind=ind+1
4546             do j=1,3
4547               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4548             enddo
4549           endif
4550         enddo
4551       enddo
4552 ! Calculate deviations from standard virtual-bond lengths
4553       ind=0
4554       do i=nnt,nct-1
4555         ind=ind+1
4556         x(ind)=vbld(i+1)**2-vbl**2
4557       enddo
4558       do i=nnt,nct
4559         if (itype(i).ne.10) then
4560           ind=ind+1
4561           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4562         endif
4563       enddo
4564       if (lprn) then
4565         write (iout,*) "Coordinates and violations"
4566         do i=1,nbond
4567           write(iout,'(i5,3f10.5,5x,e15.5)') &
4568            i,(dC_uncor(j,i),j=1,3),x(i)
4569         enddo
4570         write (iout,*) "Velocities and violations"
4571         ind=0
4572         do i=nnt,nct-1
4573           ind=ind+1
4574           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4575            i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4576         enddo
4577         do i=nnt,nct
4578           if (itype(i).ne.10) then
4579             ind=ind+1
4580             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4581              i+nres,ind,(d_t(j,i+nres),j=1,3),&
4582              scalar(d_t(1,i+nres),dC_old(1,i+nres))
4583           endif
4584         enddo
4585         write (iout,*) "gdc"
4586         do i=1,nbond
4587           write (iout,*) "i",i
4588           do j=1,nbond
4589             write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4590           enddo
4591         enddo
4592       endif
4593       xmax=dabs(x(1))
4594       do i=2,nbond
4595         if (dabs(x(i)).gt.xmax) then
4596           xmax=dabs(x(i))
4597         endif
4598       enddo
4599       if (xmax.lt.tol_rattle) then
4600         not_done=.false.
4601         goto 100
4602       endif
4603 ! Calculate the matrix of the system of equations
4604       do i=1,nbond
4605         do j=1,nbond
4606           Cmat(i,j)=0.0d0
4607           do k=1,3
4608             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4609           enddo
4610         enddo
4611       enddo
4612       if (lprn1) then
4613         write (iout,*) "Matrix Cmat"
4614         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4615       endif
4616       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4617 ! Add constraint term to positions
4618       ind=0
4619       do i=nnt,nct-1
4620         ind=ind+1
4621         do j=1,3
4622           xx=0.0d0
4623           do ii=1,nbond
4624             xx = xx+x(ii)*gdc(j,ind,ii)
4625           enddo
4626           xx=-0.5d0*xx
4627           d_t(j,i)=d_t(j,i)+xx/d_time
4628           dC(j,i)=dC(j,i)+xx
4629         enddo
4630       enddo
4631       do i=nnt,nct
4632         if (itype(i).ne.10) then
4633           ind=ind+1
4634           do j=1,3
4635             xx=0.0d0
4636             do ii=1,nbond
4637               xx = xx+x(ii)*gdc(j,ind,ii)
4638             enddo
4639             xx=-0.5d0*xx
4640             d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4641             dC(j,i+nres)=dC(j,i+nres)+xx
4642           enddo
4643         endif
4644       enddo
4645 ! Rebuild the chain using the new coordinates
4646       call chainbuild_cart
4647       if (lprn) then
4648         write (iout,*) "New coordinates, Lagrange multipliers,",&
4649         " and differences between actual and standard bond lengths"
4650         ind=0
4651         do i=nnt,nct-1
4652           ind=ind+1
4653           xx=vbld(i+1)**2-vbl**2
4654           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4655               i,(dC(j,i),j=1,3),x(ind),xx
4656         enddo
4657         do i=nnt,nct
4658           if (itype(i).ne.10) then
4659             ind=ind+1
4660             xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4661             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4662              i,(dC(j,i+nres),j=1,3),x(ind),xx
4663           endif
4664         enddo
4665         write (iout,*) "Velocities and violations"
4666         ind=0
4667         do i=nnt,nct-1
4668           ind=ind+1
4669           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4670            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4671         enddo
4672         do i=nnt,nct
4673           if (itype(i).ne.10) then
4674             ind=ind+1
4675             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4676              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4677              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4678           endif
4679         enddo
4680       endif
4681       enddo
4682   100 continue
4683       return
4684    10 write (iout,*) "Error - singularity in solving the system",&
4685        " of equations for Lagrange multipliers."
4686       stop
4687 #else
4688       write (iout,*) &
4689        "RATTLE inactive; use -DRATTLE option at compile time"
4690       stop
4691 #endif
4692       end subroutine rattle_brown
4693 !-----------------------------------------------------------------------------
4694 ! stochfric.F
4695 !-----------------------------------------------------------------------------
4696       subroutine friction_force
4697
4698       use energy_data
4699       use REMD_data
4700       use comm_syfek
4701 !      implicit real*8 (a-h,o-z)
4702 !      include 'DIMENSIONS'
4703 !      include 'COMMON.VAR'
4704 !      include 'COMMON.CHAIN'
4705 !      include 'COMMON.DERIV'
4706 !      include 'COMMON.GEO'
4707 !      include 'COMMON.LOCAL'
4708 !      include 'COMMON.INTERACT'
4709 !      include 'COMMON.MD'
4710 !#ifndef LANG0
4711 !      include 'COMMON.LANGEVIN'
4712 !#else
4713 !      include 'COMMON.LANGEVIN.lang0'
4714 !#endif
4715 !      include 'COMMON.IOUNITS'
4716 !el      real(kind=8),dimension(6*nres) :: gamvec       !(MAXRES6) maxres6=6*maxres
4717 !el      common /syfek/ gamvec
4718       real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4719 !el       ginvfric(2*nres,2*nres)       !maxres2=2*maxres
4720 !el      common /przechowalnia/ ginvfric
4721
4722       logical :: lprn = .false., checkmode = .false.
4723       integer :: i,j,ind,k,nres2,nres6
4724       nres2=2*nres
4725       nres6=6*nres
4726
4727       if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4728       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4729       do i=0,nres2
4730         do j=1,3
4731           friction(j,i)=0.0d0
4732         enddo
4733       enddo
4734
4735       do j=1,3
4736         d_t_work(j)=d_t(j,0)
4737       enddo
4738       ind=3
4739       do i=nnt,nct-1
4740         do j=1,3
4741           d_t_work(ind+j)=d_t(j,i)
4742         enddo
4743         ind=ind+3
4744       enddo
4745       do i=nnt,nct
4746         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4747           do j=1,3
4748             d_t_work(ind+j)=d_t(j,i+nres)
4749           enddo
4750           ind=ind+3
4751         endif
4752       enddo
4753
4754       call fricmat_mult(d_t_work,fric_work)
4755
4756       if (.not.checkmode) return
4757
4758       if (lprn) then
4759         write (iout,*) "d_t_work and fric_work"
4760         do i=1,3*dimen
4761           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4762         enddo
4763       endif
4764       do j=1,3
4765         friction(j,0)=fric_work(j)
4766       enddo
4767       ind=3
4768       do i=nnt,nct-1
4769         do j=1,3
4770           friction(j,i)=fric_work(ind+j)
4771         enddo
4772         ind=ind+3
4773       enddo
4774       do i=nnt,nct
4775         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4776           do j=1,3
4777             friction(j,i+nres)=fric_work(ind+j)
4778           enddo
4779           ind=ind+3
4780         endif
4781       enddo
4782       if (lprn) then
4783         write(iout,*) "Friction backbone"
4784         do i=0,nct-1
4785           write(iout,'(i5,3e15.5,5x,3e15.5)') &
4786            i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4787         enddo
4788         write(iout,*) "Friction side chain"
4789         do i=nnt,nct
4790           write(iout,'(i5,3e15.5,5x,3e15.5)') &
4791            i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4792         enddo
4793       endif
4794       if (lprn) then
4795         do j=1,3
4796           vv(j)=d_t(j,0)
4797         enddo
4798         do i=nnt,nct
4799           do j=1,3
4800             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4801             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4802             vv(j)=vv(j)+d_t(j,i)
4803           enddo
4804         enddo
4805         write (iout,*) "vvtot backbone and sidechain"
4806         do i=nnt,nct
4807           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4808            (vvtot(j,i+nres),j=1,3)
4809         enddo
4810         ind=0
4811         do i=nnt,nct-1
4812           do j=1,3
4813             v_work(ind+j)=vvtot(j,i)
4814           enddo
4815           ind=ind+3
4816         enddo
4817         do i=nnt,nct
4818           do j=1,3
4819             v_work(ind+j)=vvtot(j,i+nres)
4820           enddo
4821           ind=ind+3
4822         enddo
4823         write (iout,*) "v_work gamvec and site-based friction forces"
4824         do i=1,dimen1
4825           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4826             gamvec(i)*v_work(i)
4827         enddo
4828 !        do i=1,dimen
4829 !          fric_work1(i)=0.0d0
4830 !          do j=1,dimen1
4831 !            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
4832 !          enddo
4833 !        enddo
4834 !        write (iout,*) "fric_work and fric_work1"
4835 !        do i=1,dimen
4836 !          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
4837 !        enddo
4838         do i=1,dimen
4839           do j=1,dimen
4840             ginvfric(i,j)=0.0d0
4841             do k=1,dimen
4842               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
4843             enddo
4844           enddo
4845         enddo
4846         write (iout,*) "ginvfric"
4847         do i=1,dimen
4848           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
4849         enddo
4850         write (iout,*) "symmetry check"
4851         do i=1,dimen
4852           do j=1,i-1
4853             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
4854           enddo
4855         enddo
4856       endif
4857       return
4858       end subroutine friction_force
4859 !-----------------------------------------------------------------------------
4860       subroutine setup_fricmat
4861
4862 !     use MPI
4863       use energy_data
4864       use control_data, only:time_Bcast
4865       use control, only:tcpu
4866       use comm_syfek
4867 !      implicit real*8 (a-h,o-z)
4868 #ifdef MPI
4869       use MPI_data
4870       include 'mpif.h'
4871       real(kind=8) :: time00
4872 #endif
4873 !      include 'DIMENSIONS'
4874 !      include 'COMMON.VAR'
4875 !      include 'COMMON.CHAIN'
4876 !      include 'COMMON.DERIV'
4877 !      include 'COMMON.GEO'
4878 !      include 'COMMON.LOCAL'
4879 !      include 'COMMON.INTERACT'
4880 !      include 'COMMON.MD'
4881 !      include 'COMMON.SETUP'
4882 !      include 'COMMON.TIME1'
4883 !      integer licznik /0/
4884 !      save licznik
4885 !#ifndef LANG0
4886 !      include 'COMMON.LANGEVIN'
4887 !#else
4888 !      include 'COMMON.LANGEVIN.lang0'
4889 !#endif
4890 !      include 'COMMON.IOUNITS'
4891       integer :: IERROR
4892       integer :: i,j,ind,ind1,m
4893       logical :: lprn = .false.
4894       real(kind=8) :: dtdi !el ,gamvec(2*nres)
4895 !el      real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
4896       real(kind=8),dimension(2*nres,2*nres) :: fcopy
4897 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
4898 !el      common /syfek/ gamvec
4899       real(kind=8) :: work(8*2*nres)
4900       integer :: iwork(2*nres)
4901 !el      common /przechowalnia/ ginvfric,Ghalf,fcopy
4902       integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
4903 #ifdef MPI
4904       if (fg_rank.ne.king) goto 10
4905 #endif
4906       nres2=2*nres
4907       nres6=6*nres
4908
4909       if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
4910       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4911 !el      if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4912 !el      allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4913       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
4914
4915 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
4916 !  Zeroing out fricmat
4917       do i=1,dimen
4918         do j=1,dimen
4919           fricmat(i,j)=0.0d0
4920         enddo
4921       enddo
4922 !  Load the friction coefficients corresponding to peptide groups
4923       ind1=0
4924       do i=nnt,nct-1
4925         ind1=ind1+1
4926         gamvec(ind1)=gamp
4927       enddo
4928 !  Load the friction coefficients corresponding to side chains
4929       m=nct-nnt
4930       ind=0
4931       gamsc(ntyp1)=1.0d0
4932       do i=nnt,nct
4933         ind=ind+1
4934         ii = ind+m
4935         iti=itype(i)
4936         gamvec(ii)=gamsc(iabs(iti))
4937       enddo
4938       if (surfarea) call sdarea(gamvec)
4939 !      if (lprn) then
4940 !        write (iout,*) "Matrix A and vector gamma"
4941 !        do i=1,dimen1
4942 !          write (iout,'(i2,$)') i
4943 !          do j=1,dimen
4944 !            write (iout,'(f4.1,$)') A(i,j)
4945 !          enddo
4946 !          write (iout,'(f8.3)') gamvec(i)
4947 !        enddo
4948 !      endif
4949       if (lprn) then
4950         write (iout,*) "Vector gamvec"
4951         do i=1,dimen1
4952           write (iout,'(i5,f10.5)') i, gamvec(i)
4953         enddo
4954       endif
4955
4956 ! The friction matrix
4957       do k=1,dimen
4958        do i=1,dimen
4959          dtdi=0.0d0
4960          do j=1,dimen1
4961            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
4962          enddo
4963          fricmat(k,i)=dtdi
4964        enddo
4965       enddo
4966
4967       if (lprn) then
4968         write (iout,'(//a)') "Matrix fricmat"
4969         call matout2(dimen,dimen,nres2,nres2,fricmat)
4970       endif
4971       if (lang.eq.2 .or. lang.eq.3) then
4972 ! Mass-scale the friction matrix if non-direct integration will be performed
4973       do i=1,dimen
4974         do j=1,dimen
4975           Ginvfric(i,j)=0.0d0
4976           do k=1,dimen
4977             do l=1,dimen
4978               Ginvfric(i,j)=Ginvfric(i,j)+ &
4979                 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
4980             enddo
4981           enddo
4982         enddo
4983       enddo
4984 ! Diagonalize the friction matrix
4985       ind=0
4986       do i=1,dimen
4987         do j=1,i
4988           ind=ind+1
4989           Ghalf(ind)=Ginvfric(i,j)
4990         enddo
4991       enddo
4992       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
4993         ierr,iwork)
4994       if (lprn) then
4995         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
4996           " mass-scaled friction matrix"
4997         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
4998       endif
4999 ! Precompute matrices for tinker stochastic integrator
5000 #ifndef LANG0
5001       do i=1,dimen
5002         do j=1,dimen
5003           mt1(i,j)=0.0d0
5004           mt2(i,j)=0.0d0
5005           do k=1,dimen
5006             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5007             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5008           enddo
5009           mt3(j,i)=mt1(i,j)
5010         enddo
5011       enddo
5012 #endif
5013       else if (lang.eq.4) then
5014 ! Diagonalize the friction matrix
5015       ind=0
5016       do i=1,dimen
5017         do j=1,i
5018           ind=ind+1
5019           Ghalf(ind)=fricmat(i,j)
5020         enddo
5021       enddo
5022       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5023         ierr,iwork)
5024       if (lprn) then
5025         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5026           " friction matrix"
5027         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5028       endif
5029 ! Determine the number of zero eigenvalues of the friction matrix
5030       nzero=max0(dimen-dimen1,0)
5031 !      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5032 !        nzero=nzero+1
5033 !      enddo
5034       write (iout,*) "Number of zero eigenvalues:",nzero
5035       do i=1,dimen
5036         do j=1,dimen
5037           fricmat(i,j)=0.0d0
5038           do k=nzero+1,dimen
5039             fricmat(i,j)=fricmat(i,j) &
5040               +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5041           enddo
5042         enddo
5043       enddo
5044       if (lprn) then
5045         write (iout,'(//a)') "Generalized inverse of fricmat"
5046         call matout(dimen,dimen,nres6,nres6,fricmat)
5047       endif
5048       endif
5049 #ifdef MPI
5050   10  continue
5051       if (nfgtasks.gt.1) then
5052         if (fg_rank.eq.0) then
5053 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5054 #ifdef MPI
5055           time00=MPI_Wtime()
5056 #else
5057           time00=tcpu()
5058 #endif
5059           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5060 #ifdef MPI
5061           time_Bcast=time_Bcast+MPI_Wtime()-time00
5062 #else
5063           time_Bcast=time_Bcast+tcpu()-time00
5064 #endif
5065 !          print *,"Processor",myrank,
5066 !     &       " BROADCAST iorder in SETUP_FRICMAT"
5067         endif
5068 !       licznik=licznik+1
5069         write (iout,*) "setup_fricmat licznik"!,licznik !sp
5070 #ifdef MPI
5071         time00=MPI_Wtime()
5072 #else
5073         time00=tcpu()
5074 #endif
5075 ! Scatter the friction matrix
5076         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5077           nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5078           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5079 #ifdef TIMING
5080 #ifdef MPI
5081         time_scatter=time_scatter+MPI_Wtime()-time00
5082         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5083 #else
5084         time_scatter=time_scatter+tcpu()-time00
5085         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5086 #endif
5087 #endif
5088         do i=1,dimen
5089           do j=1,2*my_ng_count
5090             fricmat(j,i)=fcopy(i,j)
5091           enddo
5092         enddo
5093 !        write (iout,*) "My chunk of fricmat"
5094 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5095       endif
5096 #endif
5097       return
5098       end subroutine setup_fricmat
5099 !-----------------------------------------------------------------------------
5100       subroutine stochastic_force(stochforcvec)
5101
5102       use energy_data
5103       use random, only:anorm_distr
5104 !      implicit real*8 (a-h,o-z)
5105 !      include 'DIMENSIONS'
5106       use control, only: tcpu
5107       use control_data
5108 #ifdef MPI
5109       include 'mpif.h'
5110 #endif
5111 !      include 'COMMON.VAR'
5112 !      include 'COMMON.CHAIN'
5113 !      include 'COMMON.DERIV'
5114 !      include 'COMMON.GEO'
5115 !      include 'COMMON.LOCAL'
5116 !      include 'COMMON.INTERACT'
5117 !      include 'COMMON.MD'
5118 !      include 'COMMON.TIME1'
5119 !#ifndef LANG0
5120 !      include 'COMMON.LANGEVIN'
5121 !#else
5122 !      include 'COMMON.LANGEVIN.lang0'
5123 !#endif
5124 !      include 'COMMON.IOUNITS'
5125
5126       real(kind=8) :: x,sig,lowb,highb
5127       real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5128       real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5129       real(kind=8) :: time00
5130       logical :: lprn = .false.
5131       integer :: i,j,ind
5132
5133       do i=0,2*nres
5134         do j=1,3
5135           stochforc(j,i)=0.0d0
5136         enddo
5137       enddo
5138       x=0.0d0
5139
5140 #ifdef MPI
5141       time00=MPI_Wtime()
5142 #else
5143       time00=tcpu()
5144 #endif
5145 ! Compute the stochastic forces acting on bodies. Store in force.
5146       do i=nnt,nct-1
5147         sig=stdforcp(i)
5148         lowb=-5*sig
5149         highb=5*sig
5150         do j=1,3
5151           force(j,i)=anorm_distr(x,sig,lowb,highb)
5152         enddo
5153       enddo
5154       do i=nnt,nct
5155         sig2=stdforcsc(i)
5156         lowb2=-5*sig2
5157         highb2=5*sig2
5158         do j=1,3
5159           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5160         enddo
5161       enddo
5162 #ifdef MPI
5163       time_fsample=time_fsample+MPI_Wtime()-time00
5164 #else
5165       time_fsample=time_fsample+tcpu()-time00
5166 #endif
5167 ! Compute the stochastic forces acting on virtual-bond vectors.
5168       do j=1,3
5169         ff(j)=0.0d0
5170       enddo
5171       do i=nct-1,nnt,-1
5172         do j=1,3
5173           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5174         enddo
5175         do j=1,3
5176           ff(j)=ff(j)+force(j,i)
5177         enddo
5178         if (itype(i+1).ne.ntyp1) then
5179           do j=1,3
5180             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5181             ff(j)=ff(j)+force(j,i+nres+1)
5182           enddo
5183         endif
5184       enddo
5185       do j=1,3
5186         stochforc(j,0)=ff(j)+force(j,nnt+nres)
5187       enddo
5188       do i=nnt,nct
5189         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5190           do j=1,3
5191             stochforc(j,i+nres)=force(j,i+nres)
5192           enddo
5193         endif
5194       enddo
5195
5196       do j=1,3
5197         stochforcvec(j)=stochforc(j,0)
5198       enddo
5199       ind=3
5200       do i=nnt,nct-1
5201         do j=1,3
5202           stochforcvec(ind+j)=stochforc(j,i)
5203         enddo
5204         ind=ind+3
5205       enddo
5206       do i=nnt,nct
5207         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5208           do j=1,3
5209             stochforcvec(ind+j)=stochforc(j,i+nres)
5210           enddo
5211           ind=ind+3
5212         endif
5213       enddo
5214       if (lprn) then
5215         write (iout,*) "stochforcvec"
5216         do i=1,3*dimen
5217           write(iout,'(i5,e15.5)') i,stochforcvec(i)
5218         enddo
5219         write(iout,*) "Stochastic forces backbone"
5220         do i=0,nct-1
5221           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5222         enddo
5223         write(iout,*) "Stochastic forces side chain"
5224         do i=nnt,nct
5225           write(iout,'(i5,3e15.5)') &
5226             i,(stochforc(j,i+nres),j=1,3)
5227         enddo
5228       endif
5229
5230       if (lprn) then
5231
5232       ind=0
5233       do i=nnt,nct-1
5234         write (iout,*) i,ind
5235         do j=1,3
5236           forcvec(ind+j)=force(j,i)
5237         enddo
5238         ind=ind+3
5239       enddo
5240       do i=nnt,nct
5241         write (iout,*) i,ind
5242         do j=1,3
5243           forcvec(j+ind)=force(j,i+nres)
5244         enddo
5245         ind=ind+3
5246       enddo
5247
5248       write (iout,*) "forcvec"
5249       ind=0
5250       do i=nnt,nct-1
5251         do j=1,3
5252           write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5253             forcvec(ind+j)
5254         enddo
5255         ind=ind+3
5256       enddo
5257       do i=nnt,nct
5258         do j=1,3
5259           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5260            forcvec(ind+j)
5261         enddo
5262         ind=ind+3
5263       enddo
5264
5265       endif
5266
5267       return
5268       end subroutine stochastic_force
5269 !-----------------------------------------------------------------------------
5270       subroutine sdarea(gamvec)
5271 !
5272 ! Scale the friction coefficients according to solvent accessible surface areas
5273 ! Code adapted from TINKER
5274 ! AL 9/3/04
5275 !
5276       use energy_data
5277 !      implicit real*8 (a-h,o-z)
5278 !      include 'DIMENSIONS'
5279 !      include 'COMMON.CONTROL'
5280 !      include 'COMMON.VAR'
5281 !      include 'COMMON.MD'
5282 !#ifndef LANG0
5283 !      include 'COMMON.LANGEVIN'
5284 !#else
5285 !      include 'COMMON.LANGEVIN.lang0'
5286 !#endif
5287 !      include 'COMMON.CHAIN'
5288 !      include 'COMMON.DERIV'
5289 !      include 'COMMON.GEO'
5290 !      include 'COMMON.LOCAL'
5291 !      include 'COMMON.INTERACT'
5292 !      include 'COMMON.IOUNITS'
5293 !      include 'COMMON.NAMES'
5294       real(kind=8),dimension(2*nres) :: radius,gamvec   !(maxres2)
5295       real(kind=8),parameter :: twosix = 1.122462048309372981d0
5296       logical :: lprn = .false.
5297       real(kind=8) :: probe,area,ratio
5298       integer :: i,j,ind,iti
5299 !
5300 !     determine new friction coefficients every few SD steps
5301 !
5302 !     set the atomic radii to estimates of sigma values
5303 !
5304 !      print *,"Entered sdarea"
5305       probe = 0.0d0
5306
5307       do i=1,2*nres
5308         radius(i)=0.0d0
5309       enddo
5310 !  Load peptide group radii
5311       do i=nnt,nct-1
5312         radius(i)=pstok
5313       enddo
5314 !  Load side chain radii
5315       do i=nnt,nct
5316         iti=itype(i)
5317         radius(i+nres)=restok(iti)
5318       enddo
5319 !      do i=1,2*nres
5320 !        write (iout,*) "i",i," radius",radius(i)
5321 !      enddo
5322       do i = 1, 2*nres
5323          radius(i) = radius(i) / twosix
5324          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
5325       end do
5326 !
5327 !     scale atomic friction coefficients by accessible area
5328 !
5329       if (lprn) write (iout,*) &
5330         "Original gammas, surface areas, scaling factors, new gammas, ",&
5331         "std's of stochastic forces"
5332       ind=0
5333       do i=nnt,nct-1
5334         if (radius(i).gt.0.0d0) then
5335           call surfatom (i,area,radius)
5336           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5337           if (lprn) write (iout,'(i5,3f10.5,$)') &
5338             i,gamvec(ind+1),area,ratio
5339           do j=1,3
5340             ind=ind+1
5341             gamvec(ind) = ratio * gamvec(ind)
5342           enddo
5343           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
5344           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5345         endif
5346       enddo
5347       do i=nnt,nct
5348         if (radius(i+nres).gt.0.0d0) then
5349           call surfatom (i+nres,area,radius)
5350           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5351           if (lprn) write (iout,'(i5,3f10.5,$)') &
5352             i,gamvec(ind+1),area,ratio
5353           do j=1,3
5354             ind=ind+1
5355             gamvec(ind) = ratio * gamvec(ind)
5356           enddo
5357           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
5358           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5359         endif
5360       enddo
5361
5362       return
5363       end subroutine sdarea
5364 !-----------------------------------------------------------------------------
5365 ! surfatom.f
5366 !-----------------------------------------------------------------------------
5367 !
5368 !
5369 !     ###################################################
5370 !     ##  COPYRIGHT (C)  1996  by  Jay William Ponder  ##
5371 !     ##              All Rights Reserved              ##
5372 !     ###################################################
5373 !
5374 !     ################################################################
5375 !     ##                                                            ##
5376 !     ##  subroutine surfatom  --  exposed surface area of an atom  ##
5377 !     ##                                                            ##
5378 !     ################################################################
5379 !
5380 !
5381 !     "surfatom" performs an analytical computation of the surface
5382 !     area of a specified atom; a simplified version of "surface"
5383 !
5384 !     literature references:
5385 !
5386 !     T. J. Richmond, "Solvent Accessible Surface Area and
5387 !     Excluded Volume in Proteins", Journal of Molecular Biology,
5388 !     178, 63-89 (1984)
5389 !
5390 !     L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5391 !     Applied to Molecular Dynamics of Proteins in Solution",
5392 !     Protein Science, 1, 227-235 (1992)
5393 !
5394 !     variables and parameters:
5395 !
5396 !     ir       number of atom for which area is desired
5397 !     area     accessible surface area of the atom
5398 !     radius   radii of each of the individual atoms
5399 !
5400 !
5401       subroutine surfatom(ir,area,radius)
5402
5403 !      implicit real*8 (a-h,o-z)
5404 !      include 'DIMENSIONS'
5405 !      include 'sizes.i'
5406 !      include 'COMMON.GEO'
5407 !      include 'COMMON.IOUNITS'
5408 !      integer :: nres,
5409       integer :: nsup,nstart_sup
5410 !      double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5411 !      common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5412 !     & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5413 !     & dc_work(MAXRES6),nres,nres0
5414       integer,parameter :: maxarc=300
5415       integer :: i,j,k,m
5416       integer :: ii,ib,jb
5417       integer :: io,ir
5418       integer :: mi,ni,narc
5419       integer :: key(maxarc)
5420       integer :: intag(maxarc)
5421       integer :: intag1(maxarc)
5422       real(kind=8) :: area,arcsum
5423       real(kind=8) :: arclen,exang
5424       real(kind=8) :: delta,delta2
5425       real(kind=8) :: eps,rmove
5426       real(kind=8) :: xr,yr,zr
5427       real(kind=8) :: rr,rrsq
5428       real(kind=8) :: rplus,rminus
5429       real(kind=8) :: axx,axy,axz
5430       real(kind=8) :: ayx,ayy
5431       real(kind=8) :: azx,azy,azz
5432       real(kind=8) :: uxj,uyj,uzj
5433       real(kind=8) :: tx,ty,tz
5434       real(kind=8) :: txb,tyb,td
5435       real(kind=8) :: tr2,tr,txr,tyr
5436       real(kind=8) :: tk1,tk2
5437       real(kind=8) :: thec,the,t,tb
5438       real(kind=8) :: txk,tyk,tzk
5439       real(kind=8) :: t1,ti,tf,tt
5440       real(kind=8) :: txj,tyj,tzj
5441       real(kind=8) :: ccsq,cc,xysq
5442       real(kind=8) :: bsqk,bk,cosine
5443       real(kind=8) :: dsqj,gi,pix2
5444       real(kind=8) :: therk,dk,gk
5445       real(kind=8) :: risqk,rik
5446       real(kind=8) :: radius(2*nres)    !(maxatm) (maxatm=maxres2)
5447       real(kind=8) :: ri(maxarc),risq(maxarc)
5448       real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5449       real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5450       real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5451       real(kind=8) :: dsq(maxarc),bsq(maxarc)
5452       real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5453       real(kind=8) :: arci(maxarc),arcf(maxarc)
5454       real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5455       real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5456       real(kind=8) :: kent(maxarc),kout(maxarc)
5457       real(kind=8) :: ther(maxarc)
5458       logical :: moved,top
5459       logical :: omit(maxarc)
5460 !
5461 !      include 'sizes.i'
5462       maxatm = 2*nres   !maxres2 maxres2=2*maxres
5463       maxlight = 8*maxatm
5464       maxbnd = 2*maxatm
5465       maxang = 3*maxatm
5466       maxtors = 4*maxatm
5467 !
5468 !     zero out the surface area for the sphere of interest
5469 !
5470       area = 0.0d0
5471 !      write (2,*) "ir",ir," radius",radius(ir)
5472       if (radius(ir) .eq. 0.0d0)  return
5473 !
5474 !     set the overlap significance and connectivity shift
5475 !
5476       pix2 = 2.0d0 * pi
5477       delta = 1.0d-8
5478       delta2 = delta * delta
5479       eps = 1.0d-8
5480       moved = .false.
5481       rmove = 1.0d-8
5482 !
5483 !     store coordinates and radius of the sphere of interest
5484 !
5485       xr = c(1,ir)
5486       yr = c(2,ir)
5487       zr = c(3,ir)
5488       rr = radius(ir)
5489       rrsq = rr * rr
5490 !
5491 !     initialize values of some counters and summations
5492 !
5493    10 continue
5494       io = 0
5495       jb = 0
5496       ib = 0
5497       arclen = 0.0d0
5498       exang = 0.0d0
5499 !
5500 !     test each sphere to see if it overlaps the sphere of interest
5501 !
5502       do i = 1, 2*nres
5503          if (i.eq.ir .or. radius(i).eq.0.0d0)  goto 30
5504          rplus = rr + radius(i)
5505          tx = c(1,i) - xr
5506          if (abs(tx) .ge. rplus)  goto 30
5507          ty = c(2,i) - yr
5508          if (abs(ty) .ge. rplus)  goto 30
5509          tz = c(3,i) - zr
5510          if (abs(tz) .ge. rplus)  goto 30
5511 !
5512 !     check for sphere overlap by testing distance against radii
5513 !
5514          xysq = tx*tx + ty*ty
5515          if (xysq .lt. delta2) then
5516             tx = delta
5517             ty = 0.0d0
5518             xysq = delta2
5519          end if
5520          ccsq = xysq + tz*tz
5521          cc = sqrt(ccsq)
5522          if (rplus-cc .le. delta)  goto 30
5523          rminus = rr - radius(i)
5524 !
5525 !     check to see if sphere of interest is completely buried
5526 !
5527          if (cc-abs(rminus) .le. delta) then
5528             if (rminus .le. 0.0d0)  goto 170
5529             goto 30
5530          end if
5531 !
5532 !     check for too many overlaps with sphere of interest
5533 !
5534          if (io .ge. maxarc) then
5535             write (iout,20)
5536    20       format (/,' SURFATOM  --  Increase the Value of MAXARC')
5537             stop
5538          end if
5539 !
5540 !     get overlap between current sphere and sphere of interest
5541 !
5542          io = io + 1
5543          xc1(io) = tx
5544          yc1(io) = ty
5545          zc1(io) = tz
5546          dsq1(io) = xysq
5547          bsq1(io) = ccsq
5548          b1(io) = cc
5549          gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5550          intag1(io) = i
5551          omit(io) = .false.
5552    30    continue
5553       end do
5554 !
5555 !     case where no other spheres overlap the sphere of interest
5556 !
5557       if (io .eq. 0) then
5558          area = 4.0d0 * pi * rrsq
5559          return
5560       end if
5561 !
5562 !     case where only one sphere overlaps the sphere of interest
5563 !
5564       if (io .eq. 1) then
5565          area = pix2 * (1.0d0 + gr(1))
5566          area = mod(area,4.0d0*pi) * rrsq
5567          return
5568       end if
5569 !
5570 !     case where many spheres intersect the sphere of interest;
5571 !     sort the intersecting spheres by their degree of overlap
5572 !
5573       call sort2 (io,gr,key)
5574       do i = 1, io
5575          k = key(i)
5576          intag(i) = intag1(k)
5577          xc(i) = xc1(k)
5578          yc(i) = yc1(k)
5579          zc(i) = zc1(k)
5580          dsq(i) = dsq1(k)
5581          b(i) = b1(k)
5582          bsq(i) = bsq1(k)
5583       end do
5584 !
5585 !     get radius of each overlap circle on surface of the sphere
5586 !
5587       do i = 1, io
5588          gi = gr(i) * rr
5589          bg(i) = b(i) * gi
5590          risq(i) = rrsq - gi*gi
5591          ri(i) = sqrt(risq(i))
5592          ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5593       end do
5594 !
5595 !     find boundary of inaccessible area on sphere of interest
5596 !
5597       do k = 1, io-1
5598          if (.not. omit(k)) then
5599             txk = xc(k)
5600             tyk = yc(k)
5601             tzk = zc(k)
5602             bk = b(k)
5603             therk = ther(k)
5604 !
5605 !     check to see if J circle is intersecting K circle;
5606 !     get distance between circle centers and sum of radii
5607 !
5608             do j = k+1, io
5609                if (omit(j))  goto 60
5610                cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5611                cc = acos(min(1.0d0,max(-1.0d0,cc)))
5612                td = therk + ther(j)
5613 !
5614 !     check to see if circles enclose separate regions
5615 !
5616                if (cc .ge. td)  goto 60
5617 !
5618 !     check for circle J completely inside circle K
5619 !
5620                if (cc+ther(j) .lt. therk)  goto 40
5621 !
5622 !     check for circles that are essentially parallel
5623 !
5624                if (cc .gt. delta)  goto 50
5625    40          continue
5626                omit(j) = .true.
5627                goto 60
5628 !
5629 !     check to see if sphere of interest is completely buried
5630 !
5631    50          continue
5632                if (pix2-cc .le. td)  goto 170
5633    60          continue
5634             end do
5635          end if
5636       end do
5637 !
5638 !     find T value of circle intersections
5639 !
5640       do k = 1, io
5641          if (omit(k))  goto 110
5642          omit(k) = .true.
5643          narc = 0
5644          top = .false.
5645          txk = xc(k)
5646          tyk = yc(k)
5647          tzk = zc(k)
5648          dk = sqrt(dsq(k))
5649          bsqk = bsq(k)
5650          bk = b(k)
5651          gk = gr(k) * rr
5652          risqk = risq(k)
5653          rik = ri(k)
5654          therk = ther(k)
5655 !
5656 !     rotation matrix elements
5657 !
5658          t1 = tzk / (bk*dk)
5659          axx = txk * t1
5660          axy = tyk * t1
5661          axz = dk / bk
5662          ayx = tyk / dk
5663          ayy = txk / dk
5664          azx = txk / bk
5665          azy = tyk / bk
5666          azz = tzk / bk
5667          do j = 1, io
5668             if (.not. omit(j)) then
5669                txj = xc(j)
5670                tyj = yc(j)
5671                tzj = zc(j)
5672 !
5673 !     rotate spheres so K vector colinear with z-axis
5674 !
5675                uxj = txj*axx + tyj*axy - tzj*axz
5676                uyj = tyj*ayy - txj*ayx
5677                uzj = txj*azx + tyj*azy + tzj*azz
5678                cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5679                if (acos(cosine) .lt. therk+ther(j)) then
5680                   dsqj = uxj*uxj + uyj*uyj
5681                   tb = uzj*gk - bg(j)
5682                   txb = uxj * tb
5683                   tyb = uyj * tb
5684                   td = rik * dsqj
5685                   tr2 = risqk*dsqj - tb*tb
5686                   tr2 = max(eps,tr2)
5687                   tr = sqrt(tr2)
5688                   txr = uxj * tr
5689                   tyr = uyj * tr
5690 !
5691 !     get T values of intersection for K circle
5692 !
5693                   tb = (txb+tyr) / td
5694                   tb = min(1.0d0,max(-1.0d0,tb))
5695                   tk1 = acos(tb)
5696                   if (tyb-txr .lt. 0.0d0)  tk1 = pix2 - tk1
5697                   tb = (txb-tyr) / td
5698                   tb = min(1.0d0,max(-1.0d0,tb))
5699                   tk2 = acos(tb)
5700                   if (tyb+txr .lt. 0.0d0)  tk2 = pix2 - tk2
5701                   thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5702                   if (abs(thec) .lt. 1.0d0) then
5703                      the = -acos(thec)
5704                   else if (thec .ge. 1.0d0) then
5705                      the = 0.0d0
5706                   else if (thec .le. -1.0d0) then
5707                      the = -pi
5708                   end if
5709 !
5710 !     see if "tk1" is entry or exit point; check t=0 point;
5711 !     "ti" is exit point, "tf" is entry point
5712 !
5713                   cosine = min(1.0d0,max(-1.0d0, &
5714                                   (uzj*gk-uxj*rik)/(b(j)*rr)))
5715                   if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5716                      ti = tk2
5717                      tf = tk1
5718                   else
5719                      ti = tk2
5720                      tf = tk1
5721                   end if
5722                   narc = narc + 1
5723                   if (narc .ge. maxarc) then
5724                      write (iout,70)
5725    70                format (/,' SURFATOM  --  Increase the Value',&
5726                                 ' of MAXARC')
5727                      stop
5728                   end if
5729                   if (tf .le. ti) then
5730                      arcf(narc) = tf
5731                      arci(narc) = 0.0d0
5732                      tf = pix2
5733                      lt(narc) = j
5734                      ex(narc) = the
5735                      top = .true.
5736                      narc = narc + 1
5737                   end if
5738                   arcf(narc) = tf
5739                   arci(narc) = ti
5740                   lt(narc) = j
5741                   ex(narc) = the
5742                   ux(j) = uxj
5743                   uy(j) = uyj
5744                   uz(j) = uzj
5745                end if
5746             end if
5747          end do
5748          omit(k) = .false.
5749 !
5750 !     special case; K circle without intersections
5751 !
5752          if (narc .le. 0)  goto 90
5753 !
5754 !     general case; sum up arclength and set connectivity code
5755 !
5756          call sort2 (narc,arci,key)
5757          arcsum = arci(1)
5758          mi = key(1)
5759          t = arcf(mi)
5760          ni = mi
5761          if (narc .gt. 1) then
5762             do j = 2, narc
5763                m = key(j)
5764                if (t .lt. arci(j)) then
5765                   arcsum = arcsum + arci(j) - t
5766                   exang = exang + ex(ni)
5767                   jb = jb + 1
5768                   if (jb .ge. maxarc) then
5769                      write (iout,80)
5770    80                format (/,' SURFATOM  --  Increase the Value',&
5771                                 ' of MAXARC')
5772                      stop
5773                   end if
5774                   i = lt(ni)
5775                   kent(jb) = maxarc*i + k
5776                   i = lt(m)
5777                   kout(jb) = maxarc*k + i
5778                end if
5779                tt = arcf(m)
5780                if (tt .ge. t) then
5781                   t = tt
5782                   ni = m
5783                end if
5784             end do
5785          end if
5786          arcsum = arcsum + pix2 - t
5787          if (.not. top) then
5788             exang = exang + ex(ni)
5789             jb = jb + 1
5790             i = lt(ni)
5791             kent(jb) = maxarc*i + k
5792             i = lt(mi)
5793             kout(jb) = maxarc*k + i
5794          end if
5795          goto 100
5796    90    continue
5797          arcsum = pix2
5798          ib = ib + 1
5799   100    continue
5800          arclen = arclen + gr(k)*arcsum
5801   110    continue
5802       end do
5803       if (arclen .eq. 0.0d0)  goto 170
5804       if (jb .eq. 0)  goto 150
5805 !
5806 !     find number of independent boundaries and check connectivity
5807 !
5808       j = 0
5809       do k = 1, jb
5810          if (kout(k) .ne. 0) then
5811             i = k
5812   120       continue
5813             m = kout(i)
5814             kout(i) = 0
5815             j = j + 1
5816             do ii = 1, jb
5817                if (m .eq. kent(ii)) then
5818                   if (ii .eq. k) then
5819                      ib = ib + 1
5820                      if (j .eq. jb)  goto 150
5821                      goto 130
5822                   end if
5823                   i = ii
5824                   goto 120
5825                end if
5826             end do
5827   130       continue
5828          end if
5829       end do
5830       ib = ib + 1
5831 !
5832 !     attempt to fix connectivity error by moving atom slightly
5833 !
5834       if (moved) then
5835          write (iout,140)  ir
5836   140    format (/,' SURFATOM  --  Connectivity Error at Atom',i6)
5837       else
5838          moved = .true.
5839          xr = xr + rmove
5840          yr = yr + rmove
5841          zr = zr + rmove
5842          goto 10
5843       end if
5844 !
5845 !     compute the exposed surface area for the sphere of interest
5846 !
5847   150 continue
5848       area = ib*pix2 + exang + arclen
5849       area = mod(area,4.0d0*pi) * rrsq
5850 !
5851 !     attempt to fix negative area by moving atom slightly
5852 !
5853       if (area .lt. 0.0d0) then
5854          if (moved) then
5855             write (iout,160)  ir
5856   160       format (/,' SURFATOM  --  Negative Area at Atom',i6)
5857          else
5858             moved = .true.
5859             xr = xr + rmove
5860             yr = yr + rmove
5861             zr = zr + rmove
5862             goto 10
5863          end if
5864       end if
5865   170 continue
5866       return
5867       end subroutine surfatom
5868 !----------------------------------------------------------------
5869 !----------------------------------------------------------------
5870       subroutine alloc_MD_arrays
5871 !EL Allocation of arrays used by MD module
5872
5873       integer :: nres2,nres6
5874       nres2=nres*2
5875       nres6=nres*6
5876 !----------------------
5877 #ifndef LANG0
5878 ! commom.langevin
5879 !      common /langforc/
5880       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5881       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5882       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5883       allocate(fricvec(nres2,nres2))
5884       allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
5885       allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
5886       allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
5887       allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
5888       allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
5889       allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
5890       allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
5891       allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
5892       allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
5893       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5894 !      common /langmat/
5895       allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
5896 !----------------------
5897 #else
5898 ! commom.langevin.lang0
5899 !      common /langforc/
5900       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5901       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5902       allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
5903       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5904       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5905 #endif
5906
5907 !el      if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
5908 !----------------------
5909 ! commom.hairpin in CSA module
5910 !----------------------
5911 ! common.mce in MCM_MD module
5912 !----------------------
5913 ! common.MD
5914 !      common /mdgrad/ in module.energy
5915 !      common /back_constr/ in module.energy
5916 !      common /qmeas/ in module.energy
5917 !      common /mdpar/
5918 !      common /MDcalc/
5919       allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
5920 !      common /lagrange/
5921       allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
5922       allocate(d_a_work(nres6)) !(6*MAXRES)
5923 #ifdef FIVEDIAG
5924       allocate(DM(nres2),DU1(nres2),DU2(nres2))
5925       allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
5926 #else
5927       allocate(Gmat(nres2,nres2),A(nres2,nres2))
5928       if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
5929       allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
5930 #endif
5931       allocate(Geigen(nres2)) !(maxres2)
5932       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
5933 !      common /inertia/ in io_conf: parmread
5934 !      real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
5935 !      common /langevin/in io read_MDpar
5936 !      real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
5937 !      real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
5938 ! in io_conf: parmread
5939 !      real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
5940 !      common /mdpmpi/ in control: ergastulum
5941       if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
5942       if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
5943       if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
5944       if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
5945 !----------------------
5946 ! common.muca in read_muca
5947 !      common /double_muca/
5948 !      real(kind=8) :: elow,ehigh,factor,hbin,factor_min
5949 !      real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
5950 !       nemuca2,hist !(4*maxres)
5951 !      common /integer_muca/
5952 !      integer :: nmuca,imtime,muca_smooth
5953 !      common /mucarem/
5954 !      real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
5955 !----------------------
5956 ! common.MD
5957 !      common /mdgrad/ in module.energy
5958 !      common /back_constr/ in module.energy
5959 !      common /qmeas/ in module.energy
5960 !      common /mdpar/
5961 !      common /MDcalc/
5962 !      common /lagrange/
5963       allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
5964       allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
5965       allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
5966       allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
5967 !----------------------
5968 !      COMMON /BANII/ D
5969       allocate(D_ban(nres6))    !(MAXRES6) maxres6=6*maxres
5970 !      common /stochcalc/ stochforcvec
5971       allocate(stochforcvec(nres6))     !(MAXRES6) maxres6=6*maxres
5972 !----------------------
5973       return
5974       end subroutine alloc_MD_arrays
5975 !-----------------------------------------------------------------------------
5976 !-----------------------------------------------------------------------------
5977       end module MDyn