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