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