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