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.99e18 .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       integer :: iranmin,itrial,itmp
2324 !      include 'COMMON.SETUP'
2325 !      include 'COMMON.CONTROL'
2326 !      include 'COMMON.VAR'
2327 !      include 'COMMON.MD'
2328 !#ifndef LANG0
2329 !      include 'COMMON.LANGEVIN'
2330 !#else
2331 !      include 'COMMON.LANGEVIN.lang0'
2332 !#endif
2333 !      include 'COMMON.CHAIN'
2334 !      include 'COMMON.DERIV'
2335 !      include 'COMMON.GEO'
2336 !      include 'COMMON.LOCAL'
2337 !      include 'COMMON.INTERACT'
2338 !      include 'COMMON.IOUNITS'
2339 !      include 'COMMON.NAMES'
2340 !      include 'COMMON.REMD'
2341       real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2342       real(kind=8),dimension(3) :: vcm,incr,L
2343       real(kind=8) :: xv,sigv,lowb,highb
2344       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
2345       character(len=256) :: qstr
2346 !el      integer ilen
2347 !el      external ilen
2348       character(len=50) :: tytul
2349       logical :: file_exist
2350 !el      common /gucio/ cm
2351       integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2352       real(kind=8) :: etot,tt0
2353       logical :: fail
2354
2355       d_time0=d_time
2356 !      write(iout,*) "d_time", d_time
2357 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2358 ! if the friction coefficients do not depend on surface area
2359       if (lang.gt.0 .and. .not.surfarea) then
2360         do i=nnt,nct-1
2361           mnum=molnum(i)
2362           stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2363         enddo
2364         do i=nnt,nct
2365           mnum=molnum(i)
2366           stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2367                       *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2368         enddo
2369       endif
2370 ! Open the pdb file for snapshotshots
2371 #ifdef MPI
2372       if(mdpdb) then
2373         if (ilen(tmpdir).gt.0) &
2374           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2375             liczba(:ilen(liczba))//".pdb")
2376         open(ipdb,&
2377         file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2378         //".pdb")
2379       else
2380 #ifdef NOXDR
2381         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2382           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2383             liczba(:ilen(liczba))//".x")
2384         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2385         //".x"
2386 #else
2387         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2388           call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2389             liczba(:ilen(liczba))//".cx")
2390         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2391         //".cx"
2392 #endif
2393       endif
2394 #else
2395       if(mdpdb) then
2396          if (ilen(tmpdir).gt.0) &
2397            call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2398          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2399       else
2400          if (ilen(tmpdir).gt.0) &
2401            call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2402          cartname=prefix(:ilen(prefix))//"_MD.cx"
2403       endif
2404 #endif
2405       if (usampl) then
2406         write (qstr,'(256(1h ))')
2407         ipos=1
2408         do i=1,nfrag
2409           iq = qinfrag(i,iset)*10
2410           iw = wfrag(i,iset)/100
2411           if (iw.gt.0) then
2412             if(me.eq.king.or..not.out1file) &
2413              write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2414             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2415             ipos=ipos+7
2416           endif
2417         enddo
2418         do i=1,npair
2419           iq = qinpair(i,iset)*10
2420           iw = wpair(i,iset)/100
2421           if (iw.gt.0) then
2422             if(me.eq.king.or..not.out1file) &
2423              write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2424             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2425             ipos=ipos+7
2426           endif
2427         enddo
2428 !        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2429 #ifdef NOXDR
2430 !        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2431 #else
2432 !        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2433 #endif
2434 !        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2435       endif
2436       icg=1
2437       if (rest) then
2438        if (restart1file) then
2439          if (me.eq.king) &
2440            inquire(file=mremd_rst_name,exist=file_exist)
2441            write (*,*) me," Before broadcast: file_exist",file_exist
2442 #ifdef MPI !el
2443          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2444                 IERR)
2445 #endif !el
2446          write (*,*) me," After broadcast: file_exist",file_exist
2447 !        inquire(file=mremd_rst_name,exist=file_exist)
2448         if(me.eq.king.or..not.out1file) &
2449          write(iout,*) "Initial state read by master and distributed"
2450        else
2451          if (ilen(tmpdir).gt.0) &
2452            call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2453             //liczba(:ilen(liczba))//'.rst')
2454         inquire(file=rest2name,exist=file_exist)
2455        endif
2456        if(file_exist) then
2457          if(.not.restart1file) then
2458            if(me.eq.king.or..not.out1file) &
2459             write(iout,*) "Initial state will be read from file ",&
2460             rest2name(:ilen(rest2name))
2461            call readrst
2462          endif
2463          call rescale_weights(t_bath)
2464        else
2465         if(me.eq.king.or..not.out1file)then
2466          if (restart1file) then
2467           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2468              " does not exist"
2469          else
2470           write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2471              " does not exist"
2472          endif
2473          write(iout,*) "Initial velocities randomly generated"
2474         endif
2475         call random_vel
2476         totT=0.0d0
2477         totTafm=totT
2478        endif
2479       else
2480 ! Generate initial velocities
2481         if(me.eq.king.or..not.out1file) &
2482          write(iout,*) "Initial velocities randomly generated"
2483         call random_vel
2484         totT=0.0d0
2485         totTafm=totT
2486       endif
2487 !      rest2name = prefix(:ilen(prefix))//'.rst'
2488       if(me.eq.king.or..not.out1file)then
2489        write (iout,*) "Initial velocities"
2490        do i=0,nres
2491          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2492          (d_t(j,i+nres),j=1,3)
2493        enddo
2494 !  Zeroing the total angular momentum of the system
2495        write(iout,*) "Calling the zero-angular momentum subroutine"
2496       endif
2497       call inertia_tensor
2498 !  Getting the potential energy and forces and velocities and accelerations
2499       call vcm_vel(vcm)
2500 !      write (iout,*) "velocity of the center of the mass:"
2501 !      write (iout,*) (vcm(j),j=1,3)
2502       do j=1,3
2503         d_t(j,0)=d_t(j,0)-vcm(j)
2504       enddo
2505 ! Removing the velocity of the center of mass
2506       call vcm_vel(vcm)
2507       if(me.eq.king.or..not.out1file)then
2508        write (iout,*) "vcm right after adjustment:"
2509        write (iout,*) (vcm(j),j=1,3)
2510       endif
2511       if (.not.rest) then
2512 !         call chainbuild
2513          if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2514           if (overlapsc) then
2515            print *, 'Calling OVERLAP_SC'
2516            call overlap_sc(fail)
2517            print *,'after OVERLAP'
2518           endif
2519           if (searchsc) then
2520            print *,'call SC_MOVE'
2521            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2522            print *,'SC_move',nft_sc,etot
2523            if(me.eq.king.or..not.out1file) &
2524             write(iout,*) 'SC_move',nft_sc,etot
2525           endif
2526
2527           if(dccart)then
2528            print *, 'Calling MINIM_DC'
2529            call minim_dc(etot,iretcode,nfun)
2530           else
2531            call geom_to_var(nvar,varia)
2532            print *,'Calling MINIMIZE.'
2533            call minimize(etot,varia,iretcode,nfun)
2534            call var_to_geom(nvar,varia)
2535           endif
2536           if(me.eq.king.or..not.out1file) &
2537              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2538          endif
2539          if(iranconf.ne.0) then
2540 !c 8/22/17 AL Loop to produce a low-energy random conformation
2541           DO iranmin=1,40
2542           if (overlapsc) then
2543            if(me.eq.king.or..not.out1file) &
2544              write (iout,*) 'Calling OVERLAP_SC'
2545            call overlap_sc(fail)
2546           endif !endif overlap
2547
2548           if (searchsc) then
2549            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2550            print *,'SC_move',nft_sc,etot
2551            if(me.eq.king.or..not.out1file) &
2552            write(iout,*) 'SC_move',nft_sc,etot
2553           endif
2554
2555           if(dccart)then
2556            print *, 'Calling MINIM_DC'
2557            call minim_dc(etot,iretcode,nfun)
2558           else
2559            call geom_to_var(nvar,varia)
2560            print *,'Calling MINIMIZE.'
2561            call minimize(etot,varia,iretcode,nfun)
2562            call var_to_geom(nvar,varia)
2563           endif
2564           if(me.eq.king.or..not.out1file) &
2565             write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2566
2567           if (isnan(etot) .or. etot.gt.4.0d6) then
2568             write (iout,*) "Energy too large",etot, &
2569              " trying another random conformation"
2570             do itrial=1,100
2571               itmp=1
2572               call gen_rand_conf(itmp,*30)
2573               goto 40
2574    30         write (iout,*) 'Failed to generate random conformation', &
2575                ', itrial=',itrial
2576               write (*,*) 'Processor:',me, &
2577                ' Failed to generate random conformation',&
2578                ' itrial=',itrial
2579               call intout
2580 #ifdef AIX
2581               call flush_(iout)
2582 #else
2583               call flush(iout)
2584 #endif
2585             enddo
2586             write (iout,'(a,i3,a)') 'Processor:',me, &
2587              ' error in generating random conformation.'
2588             write (*,'(a,i3,a)') 'Processor:',me, &
2589              ' error in generating random conformation.'
2590             call flush(iout)
2591 #ifdef MPI
2592 !            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2593             call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2594 #else
2595             stop
2596 #endif
2597    40       continue
2598           else
2599             goto 44
2600           endif
2601           ENDDO
2602
2603           write (iout,'(a,i3,a)') 'Processor:',me, &
2604              ' failed to generate a low-energy random conformation.'
2605             write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2606              ' failed to generate a low-energy random conformation.',etot
2607             call flush(iout)
2608             call intout
2609 #ifdef MPI
2610 !            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2611         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2612 #else
2613             stop
2614 #endif
2615    44     continue
2616         endif
2617       endif
2618       call chainbuild_cart
2619       call kinetic(EK)
2620       if (tbf) then
2621         call verlet_bath
2622       endif
2623       kinetic_T=2.0d0/(dimen3*Rb)*EK
2624       if(me.eq.king.or..not.out1file)then
2625        call cartprint
2626        call intout
2627       endif
2628 #ifdef MPI
2629       tt0=MPI_Wtime()
2630 #else
2631       tt0=tcpu()
2632 #endif
2633       call zerograd
2634       write(iout,*) "before ETOTAL"
2635       call etotal(potEcomp)
2636       if (large) call enerprint(potEcomp)
2637 #ifdef TIMING_ENE
2638 #ifdef MPI
2639       t_etotal=t_etotal+MPI_Wtime()-tt0
2640 #else
2641       t_etotal=t_etotal+tcpu()-tt0
2642 #endif
2643 #endif
2644       potE=potEcomp(0)
2645       call cartgrad
2646       call lagrangian
2647       call max_accel
2648       if (amax*d_time .gt. dvmax) then
2649         d_time=d_time*dvmax/amax
2650         if(me.eq.king.or..not.out1file) write (iout,*) &
2651          "Time step reduced to",d_time,&
2652          " because of too large initial acceleration."
2653       endif
2654       if(me.eq.king.or..not.out1file)then
2655        write(iout,*) "Potential energy and its components"
2656        call enerprint(potEcomp)
2657 !       write(iout,*) (potEcomp(i),i=0,n_ene)
2658       endif
2659       potE=potEcomp(0)-potEcomp(20)
2660       totE=EK+potE
2661       itime=0
2662       if (ntwe.ne.0) call statout(itime)
2663       if(me.eq.king.or..not.out1file) &
2664         write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2665          " Kinetic energy",EK," Potential energy",potE, &
2666          " Total energy",totE," Maximum acceleration ", &
2667          amax
2668       if (large) then
2669         write (iout,*) "Initial coordinates"
2670         do i=1,nres
2671           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2672           (c(j,i+nres),j=1,3)
2673         enddo
2674         write (iout,*) "Initial dC"
2675         do i=0,nres
2676           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2677           (dc(j,i+nres),j=1,3)
2678         enddo
2679         write (iout,*) "Initial velocities"
2680         write (iout,"(13x,' backbone ',23x,' side chain')")
2681         do i=0,nres
2682           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2683           (d_t(j,i+nres),j=1,3)
2684         enddo
2685         write (iout,*) "Initial accelerations"
2686         do i=0,nres
2687 !          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2688           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2689           (d_a(j,i+nres),j=1,3)
2690         enddo
2691       endif
2692       do i=0,2*nres
2693         do j=1,3
2694           dc_old(j,i)=dc(j,i)
2695           d_t_old(j,i)=d_t(j,i)
2696           d_a_old(j,i)=d_a(j,i)
2697         enddo
2698 !        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2699       enddo
2700       if (RESPA) then
2701 #ifdef MPI
2702         tt0 =MPI_Wtime()
2703 #else
2704         tt0 = tcpu()
2705 #endif
2706         call zerograd
2707         call etotal_short(energia_short)
2708         if (large) call enerprint(potEcomp)
2709 #ifdef TIMING_ENE
2710 #ifdef MPI
2711         t_eshort=t_eshort+MPI_Wtime()-tt0
2712 #else
2713         t_eshort=t_eshort+tcpu()-tt0
2714 #endif
2715 #endif
2716         call cartgrad
2717         call lagrangian
2718         if(.not.out1file .and. large) then
2719           write (iout,*) "energia_long",energia_long(0),&
2720            " energia_short",energia_short(0),&
2721            " total",energia_long(0)+energia_short(0)
2722           write (iout,*) "Initial fast-force accelerations"
2723           do i=0,nres
2724             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2725             (d_a(j,i+nres),j=1,3)
2726           enddo
2727         endif
2728 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2729         do i=0,2*nres
2730           do j=1,3
2731             d_a_short(j,i)=d_a(j,i)
2732           enddo
2733         enddo
2734 #ifdef MPI
2735         tt0=MPI_Wtime()
2736 #else
2737         tt0=tcpu()
2738 #endif
2739         call zerograd
2740         call etotal_long(energia_long)
2741         if (large) call enerprint(potEcomp)
2742 #ifdef TIMING_ENE
2743 #ifdef MPI
2744         t_elong=t_elong+MPI_Wtime()-tt0
2745 #else
2746         t_elong=t_elong+tcpu()-tt0
2747 #endif
2748 #endif
2749         call cartgrad
2750         call lagrangian
2751         if(.not.out1file .and. large) then
2752           write (iout,*) "energia_long",energia_long(0)
2753           write (iout,*) "Initial slow-force accelerations"
2754           do i=0,nres
2755             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2756             (d_a(j,i+nres),j=1,3)
2757           enddo
2758         endif
2759 #ifdef MPI
2760         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2761 #else
2762         t_enegrad=t_enegrad+tcpu()-tt0
2763 #endif
2764       endif
2765       return
2766       end subroutine init_MD
2767 !-----------------------------------------------------------------------------
2768       subroutine random_vel
2769
2770 !      implicit real*8 (a-h,o-z)
2771       use energy_data
2772       use random, only:anorm_distr
2773       use MD_data
2774 !      include 'DIMENSIONS'
2775 !      include 'COMMON.CONTROL'
2776 !      include 'COMMON.VAR'
2777 !      include 'COMMON.MD'
2778 !#ifndef LANG0
2779 !      include 'COMMON.LANGEVIN'
2780 !#else
2781 !      include 'COMMON.LANGEVIN.lang0'
2782 !#endif
2783 !      include 'COMMON.CHAIN'
2784 !      include 'COMMON.DERIV'
2785 !      include 'COMMON.GEO'
2786 !      include 'COMMON.LOCAL'
2787 !      include 'COMMON.INTERACT'
2788 !      include 'COMMON.IOUNITS'
2789 !      include 'COMMON.NAMES'
2790 !      include 'COMMON.TIME1'
2791       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
2792 !#define DEBUG
2793 #ifdef FIVEDIAG
2794        real(kind=8) ,allocatable, dimension(:)  :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2795        real(kind=8) :: sumx
2796 #ifdef DEBUG
2797        real(kind=8) ,allocatable, dimension(:)  :: rsold
2798        real (kind=8),allocatable,dimension(:,:) :: matold
2799        integer :: iti
2800 #endif
2801 #endif
2802       integer :: i,j,ii,k,ind,mark,imark,mnum
2803 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2804 ! First generate velocities in the eigenspace of the G matrix
2805 !      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2806 !      call flush(iout)
2807       xv=0.0d0
2808       ii=0
2809       do i=1,dimen
2810         do k=1,3
2811           ii=ii+1
2812           sigv=dsqrt((Rb*t_bath)/geigen(i))
2813           lowb=-5*sigv
2814           highb=5*sigv
2815           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2816 #ifdef DEBUG
2817           write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2818             " d_t_work_new",d_t_work_new(ii)
2819 #endif
2820         enddo
2821       enddo
2822 #ifdef DEBUG
2823 ! diagnostics
2824       Ek1=0.0d0
2825       ii=0
2826       do i=1,dimen
2827         do k=1,3
2828           ii=ii+1
2829           Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2830         enddo
2831       enddo
2832       write (iout,*) "Ek from eigenvectors",Ek1
2833       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2834 ! end diagnostics
2835 #endif
2836 #ifdef FIVEDIAG
2837 ! Transform velocities to UNRES coordinate space
2838        allocate (DL1(2*nres))
2839        allocate (DDU1(2*nres))
2840        allocate (DL2(2*nres))
2841        allocate (DDU2(2*nres))
2842        allocate (xsolv(2*nres))
2843        allocate (DML(2*nres))
2844        allocate (rs(2*nres))
2845 #ifdef DEBUG
2846        allocate (rsold(2*nres))
2847        allocate (matold(2*nres,2*nres))
2848        matold=0
2849        matold(1,1)=DMorig(1)
2850        matold(1,2)=DU1orig(1)
2851        matold(1,3)=DU2orig(1)
2852        write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2853
2854       do i=2,dimen-1
2855         do j=1,dimen
2856           if (i.eq.j) then
2857              matold(i,j)=DMorig(i)
2858              matold(i,j-1)=DU1orig(i-1)
2859            if (j.ge.3) then
2860                  matold(i,j-2)=DU2orig(i-2)
2861                endif
2862
2863             endif
2864
2865           enddo
2866           do j=1,dimen-1
2867             if (i.eq.j) then
2868               matold(i,j+1)=DU1orig(i)
2869
2870              end if
2871           enddo
2872           do j=1,dimen-2
2873             if(i.eq.j) then
2874                matold(i,j+2)=DU2orig(i)
2875             endif
2876           enddo
2877        enddo
2878        matold(dimen,dimen)=DMorig(dimen)
2879        matold(dimen,dimen-1)=DU1orig(dimen-1)
2880        matold(dimen,dimen-2)=DU2orig(dimen-2)
2881        write(iout,*) "old gmatrix"
2882        call matout(dimen,dimen,2*nres,2*nres,matold)
2883 #endif
2884       d_t_work=0.0d0
2885       do i=1,dimen
2886 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2887        do imark=dimen,1,-1
2888
2889          do j=1,imark-1
2890            DML(j)=DMorig(j)-geigen(i)
2891          enddo
2892          do j=imark+1,dimen
2893            DML(j-1)=DMorig(j)-geigen(i)
2894          enddo
2895          do j=1,imark-2
2896            DDU1(j)=DU1orig(j)
2897          enddo
2898          DDU1(imark-1)=DU2orig(imark-1)
2899          do j=imark+1,dimen-1
2900            DDU1(j-1)=DU1orig(j)
2901          enddo
2902          do j=1,imark-3
2903            DDU2(j)=DU2orig(j)
2904          enddo
2905          DDU2(imark-2)=0.0d0
2906          DDU2(imark-1)=0.0d0
2907          do j=imark,dimen-3
2908            DDU2(j)=DU2orig(j+1)
2909          enddo
2910          do j=1,dimen-3
2911            DL2(j+2)=DDU2(j)
2912          enddo
2913          do j=1,dimen-2
2914            DL1(j+1)=DDU1(j)
2915          enddo
2916 #ifdef DEBUG
2917          write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2918          write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2919          write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2920          write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2921          write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2922          write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2923 #endif
2924          rs=0.0d0
2925          if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2926          if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2927          if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2928          if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2929 #ifdef DEBUG
2930          rsold=rs
2931 #endif
2932 !         write (iout,*) "Vector rs"
2933 !         do j=1,dimen-1
2934 !           write (iout,*) j,rs(j)
2935 !         enddo
2936          xsolv=0.0d0
2937          call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2938
2939          if (mark.eq.1) then
2940
2941 #ifdef DEBUG
2942 ! Check solution
2943          do j=1,imark-1
2944            sumx=-geigen(i)*xsolv(j)
2945            do k=1,imark-1
2946              sumx=sumx+matold(j,k)*xsolv(k)
2947            enddo
2948            do k=imark+1,dimen
2949              sumx=sumx+matold(j,k)*xsolv(k-1)
2950            enddo
2951            write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2952          enddo
2953          do j=imark+1,dimen
2954            sumx=-geigen(i)*xsolv(j-1)
2955            do k=1,imark-1
2956              sumx=sumx+matold(j,k)*xsolv(k)
2957            enddo
2958            do k=imark+1,dimen
2959              sumx=sumx+matold(j,k)*xsolv(k-1)
2960            enddo
2961            write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2962          enddo
2963 ! end check
2964          write (iout,*)&
2965           "Solution of equations system",i," for eigenvalue",geigen(i)
2966          do j=1,dimen-1
2967            write(iout,'(i5,f10.5)') j,xsolv(j)
2968          enddo
2969 #endif
2970          do j=dimen-1,imark,-1
2971            xsolv(j+1)=xsolv(j)
2972          enddo
2973          xsolv(imark)=1.0d0
2974 #ifdef DEBUG
2975          write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2976          do j=1,dimen
2977            write(iout,'(i5,f10.5)') j,xsolv(j)
2978          enddo
2979 #endif
2980 ! Normalize ith eigenvector
2981          sumx=0.0d0
2982          do j=1,dimen
2983            sumx=sumx+xsolv(j)**2
2984          enddo
2985          sumx=dsqrt(sumx)
2986          do j=1,dimen
2987            xsolv(j)=xsolv(j)/sumx
2988          enddo
2989 #ifdef DEBUG
2990          write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2991          do j=1,dimen
2992            write(iout,'(i5,3f10.5)') j,xsolv(j)
2993          enddo
2994 #endif
2995 ! All done at this point for eigenvector i, exit loop
2996          exit
2997
2998          endif ! mark.eq.1
2999
3000        enddo ! imark
3001
3002        if (mark.ne.1) then
3003          write (iout,*) "Unable to find eigenvector",i
3004        endif
3005
3006 !       write (iout,*) "i=",i
3007        do k=1,3
3008 !         write (iout,*) "k=",k
3009          do j=1,dimen
3010            ind=(j-1)*3+k
3011 !           write(iout,*) "ind",ind," ind1",3*(i-1)+k
3012            d_t_work(ind)=d_t_work(ind) &
3013                             +xsolv(j)*d_t_work_new(3*(i-1)+k)
3014          enddo
3015        enddo
3016       enddo ! i (loop over the eigenvectors)
3017
3018 #ifdef DEBUG
3019       write (iout,*) "d_t_work"
3020       do i=1,3*dimen
3021         write (iout,"(i5,f10.5)") i,d_t_work(i)
3022       enddo
3023       Ek1=0.0d0
3024       ii=0
3025       do i=nnt,nct
3026 !        if (itype(i,1).eq.10) then
3027          mnum=molnum(i)
3028           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3029         iti=iabs(itype(i,mnum))
3030 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3031          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3032           .or.(mnum.eq.5)) then
3033           j=ii+3
3034         else
3035           j=ii+6
3036         endif
3037         if (i.lt.nct) then
3038           do k=1,3
3039 !            write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3040             Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3041             0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3042           enddo
3043         endif
3044          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3045           .and.(mnum.ne.5)) ii=ii+3
3046         write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3047         write (iout,*) "ii",ii
3048         do k=1,3
3049           ii=ii+1
3050           write (iout,*) "k",k," ii",ii,"EK1",EK1
3051           if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3052           .and.(mnum.ne.5))&
3053            Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3054           Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3055         enddo
3056         write (iout,*) "i",i," ii",ii
3057       enddo
3058       write (iout,*) "Ek from d_t_work",Ek1
3059       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3060 #endif
3061       if(allocated(DDU1)) deallocate(DDU1)
3062       if(allocated(DDU2)) deallocate(DDU2)
3063       if(allocated(DL2)) deallocate(DL2)
3064       if(allocated(DL1)) deallocate(DL1)
3065       if(allocated(xsolv)) deallocate(xsolv)
3066       if(allocated(DML)) deallocate(DML)
3067       if(allocated(rs)) deallocate(rs)
3068 #ifdef DEBUG
3069       if(allocated(matold)) deallocate(matold)
3070       if(allocated(rsold)) deallocate(rsold)
3071 #endif
3072       ind=1
3073       do j=nnt,nct
3074         do k=1,3
3075           d_t(k,j)=d_t_work(ind)
3076           ind=ind+1
3077         enddo
3078          mnum=molnum(j)
3079          if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3080           .and.(mnum.ne.5)) then
3081           do k=1,3
3082             d_t(k,j+nres)=d_t_work(ind)
3083             ind=ind+1
3084           enddo
3085         end if
3086       enddo
3087 #ifdef DEBUG
3088       write (iout,*) "Random velocities in the Calpha,SC space"
3089       do i=nnt,nct
3090         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3091       enddo
3092       do i=nnt,nct
3093         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3094       enddo
3095 #endif
3096       do j=1,3
3097         d_t(j,0)=d_t(j,nnt)
3098       enddo
3099       do i=nnt,nct
3100 !        if (itype(i,1).eq.10) then
3101          mnum=molnum(i)
3102          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3103           .or.(mnum.eq.5)) then
3104           do j=1,3
3105             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3106           enddo
3107         else
3108           do j=1,3
3109             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3110             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3111           enddo
3112         end if
3113       enddo
3114 #ifdef DEBUG
3115       write (iout,*)"Random velocities in the virtual-bond-vector space"
3116       do i=nnt,nct-1
3117         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3118       enddo
3119       do i=nnt,nct
3120         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3121       enddo
3122       call kinetic(Ek1)
3123       write (iout,*) "Ek from d_t_work",Ek1
3124       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3125 #endif
3126 #else
3127       do k=0,2
3128         do i=1,dimen
3129           ind=(i-1)*3+k+1
3130           d_t_work(ind)=0.0d0
3131           do j=1,dimen
3132             d_t_work(ind)=d_t_work(ind) &
3133                             +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3134           enddo
3135 !          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3136 !          call flush(iout)
3137         enddo
3138       enddo
3139 ! Transfer to the d_t vector
3140       do j=1,3
3141         d_t(j,0)=d_t_work(j)
3142       enddo
3143       ind=3
3144       do i=nnt,nct-1
3145         do j=1,3
3146           ind=ind+1
3147           d_t(j,i)=d_t_work(ind)
3148         enddo
3149       enddo
3150       do i=nnt,nct
3151          mnum=molnum(i)
3152 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3153          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3154           .and.(mnum.ne.5)) then
3155           do j=1,3
3156             ind=ind+1
3157             d_t(j,i+nres)=d_t_work(ind)
3158           enddo
3159         endif
3160       enddo
3161 #endif
3162 !      call kinetic(EK)
3163 !      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3164 !        2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3165 !      call flush(iout)
3166 !      write(iout,*) "end init MD"
3167       return
3168       end subroutine random_vel
3169 !-----------------------------------------------------------------------------
3170 #ifndef LANG0
3171       subroutine sd_verlet_p_setup
3172 ! Sets up the parameters of stochastic Verlet algorithm
3173 !      implicit real*8 (a-h,o-z)
3174 !      include 'DIMENSIONS'
3175       use control, only: tcpu
3176       use control_data
3177 #ifdef MPI
3178       include 'mpif.h'
3179 #endif
3180 !      include 'COMMON.CONTROL'
3181 !      include 'COMMON.VAR'
3182 !      include 'COMMON.MD'
3183 !#ifndef LANG0
3184 !      include 'COMMON.LANGEVIN'
3185 !#else
3186 !      include 'COMMON.LANGEVIN.lang0'
3187 !#endif
3188 !      include 'COMMON.CHAIN'
3189 !      include 'COMMON.DERIV'
3190 !      include 'COMMON.GEO'
3191 !      include 'COMMON.LOCAL'
3192 !      include 'COMMON.INTERACT'
3193 !      include 'COMMON.IOUNITS'
3194 !      include 'COMMON.NAMES'
3195 !      include 'COMMON.TIME1'
3196       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
3197       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3198       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3199        prand_vec,vrand_vec1,vrand_vec2  !(MAXRES6) maxres6=6*maxres
3200       logical :: lprn = .false.
3201       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3202       real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3203                  gdt9,psig,tt0
3204       integer :: i,maxres2
3205 #ifdef MPI
3206       tt0 = MPI_Wtime()
3207 #else
3208       tt0 = tcpu()
3209 #endif
3210 !
3211 ! AL 8/17/04 Code adapted from tinker
3212 !
3213 ! Get the frictional and random terms for stochastic dynamics in the
3214 ! eigenspace of mass-scaled UNRES friction matrix
3215 !
3216       maxres2=2*nres
3217       do i = 1, dimen
3218             gdt = fricgam(i) * d_time
3219 !
3220 ! Stochastic dynamics reduces to simple MD for zero friction
3221 !
3222             if (gdt .le. zero) then
3223                pfric_vec(i) = 1.0d0
3224                vfric_vec(i) = d_time
3225                afric_vec(i) = 0.5d0 * d_time * d_time
3226                prand_vec(i) = 0.0d0
3227                vrand_vec1(i) = 0.0d0
3228                vrand_vec2(i) = 0.0d0
3229 !
3230 ! Analytical expressions when friction coefficient is large
3231 !
3232             else
3233                if (gdt .ge. gdt_radius) then
3234                   egdt = dexp(-gdt)
3235                   pfric_vec(i) = egdt
3236                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3237                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3238                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3239                   vterm = 1.0d0 - egdt**2
3240                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3241 !
3242 ! Use series expansions when friction coefficient is small
3243 !
3244                else
3245                   gdt2 = gdt * gdt
3246                   gdt3 = gdt * gdt2
3247                   gdt4 = gdt2 * gdt2
3248                   gdt5 = gdt2 * gdt3
3249                   gdt6 = gdt3 * gdt3
3250                   gdt7 = gdt3 * gdt4
3251                   gdt8 = gdt4 * gdt4
3252                   gdt9 = gdt4 * gdt5
3253                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3254                                 - gdt5/120.0d0 + gdt6/720.0d0 &
3255                                 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3256                                 - gdt9/362880.0d0) / fricgam(i)**2
3257                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3258                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3259                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3260                              + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3261                              + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3262                              + 127.0d0*gdt9/90720.0d0
3263                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3264                              - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3265                              - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3266                              - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3267                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3268                              - 17.0d0*gdt2/1280.0d0 &
3269                              + 17.0d0*gdt3/6144.0d0 &
3270                              + 40967.0d0*gdt4/34406400.0d0 &
3271                              - 57203.0d0*gdt5/275251200.0d0 &
3272                              - 1429487.0d0*gdt6/13212057600.0d0)
3273                end if
3274 !
3275 ! Compute the scaling factors of random terms for the nonzero friction case
3276 !
3277                ktm = 0.5d0*d_time/fricgam(i)
3278                psig = dsqrt(ktm*pterm) / fricgam(i)
3279                vsig = dsqrt(ktm*vterm)
3280                rhoc = dsqrt(1.0d0 - rho*rho)
3281                prand_vec(i) = psig
3282                vrand_vec1(i) = vsig * rho
3283                vrand_vec2(i) = vsig * rhoc
3284             end if
3285       end do
3286       if (lprn) then
3287       write (iout,*) &
3288         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3289         " vrand_vec2"
3290       do i=1,dimen
3291         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3292             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3293       enddo
3294       endif
3295 !
3296 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3297 !
3298 #ifndef   LANG0
3299       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3300       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3301       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3302       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3303       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3304       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3305 #endif
3306 #ifdef MPI
3307       t_sdsetup=t_sdsetup+MPI_Wtime()
3308 #else
3309       t_sdsetup=t_sdsetup+tcpu()-tt0
3310 #endif
3311       return
3312       end subroutine sd_verlet_p_setup
3313 !-----------------------------------------------------------------------------
3314       subroutine eigtransf1(n,ndim,ab,d,c)
3315
3316 !el      implicit none
3317       integer :: n,ndim
3318       real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3319       integer :: i,j,k
3320       do i=1,n
3321         do j=1,n
3322           c(i,j)=0.0d0
3323           do k=1,n
3324             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3325           enddo
3326         enddo
3327       enddo
3328       return
3329       end subroutine eigtransf1
3330 !-----------------------------------------------------------------------------
3331       subroutine eigtransf(n,ndim,a,b,d,c)
3332
3333 !el      implicit none
3334       integer :: n,ndim
3335       real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3336       integer :: i,j,k
3337       do i=1,n
3338         do j=1,n
3339           c(i,j)=0.0d0
3340           do k=1,n
3341             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3342           enddo
3343         enddo
3344       enddo
3345       return
3346       end subroutine eigtransf
3347 !-----------------------------------------------------------------------------
3348       subroutine sd_verlet1
3349
3350 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3351       use energy_data
3352 !      implicit real*8 (a-h,o-z)
3353 !      include 'DIMENSIONS'
3354 !      include 'COMMON.CONTROL'
3355 !      include 'COMMON.VAR'
3356 !      include 'COMMON.MD'
3357 !#ifndef LANG0
3358 !      include 'COMMON.LANGEVIN'
3359 !#else
3360 !      include 'COMMON.LANGEVIN.lang0'
3361 !#endif
3362 !      include 'COMMON.CHAIN'
3363 !      include 'COMMON.DERIV'
3364 !      include 'COMMON.GEO'
3365 !      include 'COMMON.LOCAL'
3366 !      include 'COMMON.INTERACT'
3367 !      include 'COMMON.IOUNITS'
3368 !      include 'COMMON.NAMES'
3369 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3370 !el      common /stochcalc/ stochforcvec
3371       logical :: lprn = .false.
3372       real(kind=8) :: ddt1,ddt2
3373       integer :: i,j,ind,inres
3374
3375 !      write (iout,*) "dc_old"
3376 !      do i=0,nres
3377 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3378 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3379 !      enddo
3380       do j=1,3
3381         dc_work(j)=dc_old(j,0)
3382         d_t_work(j)=d_t_old(j,0)
3383         d_a_work(j)=d_a_old(j,0)
3384       enddo
3385       ind=3
3386       do i=nnt,nct-1
3387         do j=1,3
3388           dc_work(ind+j)=dc_old(j,i)
3389           d_t_work(ind+j)=d_t_old(j,i)
3390           d_a_work(ind+j)=d_a_old(j,i)
3391         enddo
3392         ind=ind+3
3393       enddo
3394       do i=nnt,nct
3395          mnum=molnum(i)
3396 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3397          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3398           .and.(mnum.ne.5)) then
3399           do j=1,3
3400             dc_work(ind+j)=dc_old(j,i+nres)
3401             d_t_work(ind+j)=d_t_old(j,i+nres)
3402             d_a_work(ind+j)=d_a_old(j,i+nres)
3403           enddo
3404           ind=ind+3
3405         endif
3406       enddo
3407 #ifndef LANG0
3408       if (lprn) then
3409       write (iout,*) &
3410         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3411         " vrand_mat2"
3412       do i=1,dimen
3413         do j=1,dimen
3414           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3415             vfric_mat(i,j),afric_mat(i,j),&
3416             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3417         enddo
3418       enddo
3419       endif
3420       do i=1,dimen
3421         ddt1=0.0d0
3422         ddt2=0.0d0
3423         do j=1,dimen
3424           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3425             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3426           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3427           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3428         enddo
3429         d_t_work_new(i)=ddt1+0.5d0*ddt2
3430         d_t_work(i)=ddt1+ddt2
3431       enddo
3432 #endif
3433       do j=1,3
3434         dc(j,0)=dc_work(j)
3435         d_t(j,0)=d_t_work(j)
3436       enddo
3437       ind=3
3438       do i=nnt,nct-1
3439         do j=1,3
3440           dc(j,i)=dc_work(ind+j)
3441           d_t(j,i)=d_t_work(ind+j)
3442         enddo
3443         ind=ind+3
3444       enddo
3445       do i=nnt,nct
3446          mnum=molnum(i)
3447 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3448          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3449           .and.(mnum.ne.5)) then
3450           inres=i+nres
3451           do j=1,3
3452             dc(j,inres)=dc_work(ind+j)
3453             d_t(j,inres)=d_t_work(ind+j)
3454           enddo
3455           ind=ind+3
3456         endif
3457       enddo
3458       return
3459       end subroutine sd_verlet1
3460 !-----------------------------------------------------------------------------
3461       subroutine sd_verlet2
3462
3463 !  Calculating the adjusted velocities for accelerations
3464       use energy_data
3465 !      implicit real*8 (a-h,o-z)
3466 !      include 'DIMENSIONS'
3467 !      include 'COMMON.CONTROL'
3468 !      include 'COMMON.VAR'
3469 !      include 'COMMON.MD'
3470 !#ifndef LANG0
3471 !      include 'COMMON.LANGEVIN'
3472 !#else
3473 !      include 'COMMON.LANGEVIN.lang0'
3474 !#endif
3475 !      include 'COMMON.CHAIN'
3476 !      include 'COMMON.DERIV'
3477 !      include 'COMMON.GEO'
3478 !      include 'COMMON.LOCAL'
3479 !      include 'COMMON.INTERACT'
3480 !      include 'COMMON.IOUNITS'
3481 !      include 'COMMON.NAMES'
3482 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3483        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3484 !el      common /stochcalc/ stochforcvec
3485 !
3486       real(kind=8) :: ddt1,ddt2
3487       integer :: i,j,ind,inres
3488 ! Compute the stochastic forces which contribute to velocity change
3489 !
3490       call stochastic_force(stochforcvecV)
3491
3492 #ifndef LANG0
3493       do i=1,dimen
3494         ddt1=0.0d0
3495         ddt2=0.0d0
3496         do j=1,dimen
3497           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3498           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3499            vrand_mat2(i,j)*stochforcvecV(j)
3500         enddo
3501         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3502       enddo
3503 #endif
3504       do j=1,3
3505         d_t(j,0)=d_t_work(j)
3506       enddo
3507       ind=3
3508       do i=nnt,nct-1
3509         do j=1,3
3510           d_t(j,i)=d_t_work(ind+j)
3511         enddo
3512         ind=ind+3
3513       enddo
3514       do i=nnt,nct
3515          mnum=molnum(i)
3516 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3517          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3518           .and.(mnum.ne.5)) then
3519           inres=i+nres
3520           do j=1,3
3521             d_t(j,inres)=d_t_work(ind+j)
3522           enddo
3523           ind=ind+3
3524         endif
3525       enddo
3526       return
3527       end subroutine sd_verlet2
3528 !-----------------------------------------------------------------------------
3529       subroutine sd_verlet_ciccotti_setup
3530
3531 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3532 ! version
3533 !      implicit real*8 (a-h,o-z)
3534 !      include 'DIMENSIONS'
3535       use control, only: tcpu
3536       use control_data
3537 #ifdef MPI
3538       include 'mpif.h'
3539 #endif
3540 !      include 'COMMON.CONTROL'
3541 !      include 'COMMON.VAR'
3542 !      include 'COMMON.MD'
3543 !#ifndef LANG0
3544 !      include 'COMMON.LANGEVIN'
3545 !#else
3546 !      include 'COMMON.LANGEVIN.lang0'
3547 !#endif
3548 !      include 'COMMON.CHAIN'
3549 !      include 'COMMON.DERIV'
3550 !      include 'COMMON.GEO'
3551 !      include 'COMMON.LOCAL'
3552 !      include 'COMMON.INTERACT'
3553 !      include 'COMMON.IOUNITS'
3554 !      include 'COMMON.NAMES'
3555 !      include 'COMMON.TIME1'
3556       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
3557       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3558       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3559         prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3560       logical :: lprn = .false.
3561       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3562       real(kind=8) :: ktm,gdt,egdt,tt0
3563       integer :: i,maxres2
3564 #ifdef MPI
3565       tt0 = MPI_Wtime()
3566 #else
3567       tt0 = tcpu()
3568 #endif
3569 !
3570 ! AL 8/17/04 Code adapted from tinker
3571 !
3572 ! Get the frictional and random terms for stochastic dynamics in the
3573 ! eigenspace of mass-scaled UNRES friction matrix
3574 !
3575       maxres2=2*nres
3576       do i = 1, dimen
3577             write (iout,*) "i",i," fricgam",fricgam(i)
3578             gdt = fricgam(i) * d_time
3579 !
3580 ! Stochastic dynamics reduces to simple MD for zero friction
3581 !
3582             if (gdt .le. zero) then
3583                pfric_vec(i) = 1.0d0
3584                vfric_vec(i) = d_time
3585                afric_vec(i) = 0.5d0*d_time*d_time
3586                prand_vec(i) = afric_vec(i)
3587                vrand_vec2(i) = vfric_vec(i)
3588 !
3589 ! Analytical expressions when friction coefficient is large
3590 !
3591             else
3592                egdt = dexp(-gdt)
3593                pfric_vec(i) = egdt
3594                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3595                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3596                prand_vec(i) = afric_vec(i)
3597                vrand_vec2(i) = vfric_vec(i)
3598 !
3599 ! Compute the scaling factors of random terms for the nonzero friction case
3600 !
3601 !               ktm = 0.5d0*d_time/fricgam(i)
3602 !               psig = dsqrt(ktm*pterm) / fricgam(i)
3603 !               vsig = dsqrt(ktm*vterm)
3604 !               prand_vec(i) = psig*afric_vec(i)
3605 !               vrand_vec2(i) = vsig*vfric_vec(i)
3606             end if
3607       end do
3608       if (lprn) then
3609       write (iout,*) &
3610         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3611         " vrand_vec2"
3612       do i=1,dimen
3613         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3614             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3615       enddo
3616       endif
3617 !
3618 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3619 !
3620       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3621       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3622       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3623       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3624       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3625 #ifdef MPI
3626       t_sdsetup=t_sdsetup+MPI_Wtime()
3627 #else
3628       t_sdsetup=t_sdsetup+tcpu()-tt0
3629 #endif
3630       return
3631       end subroutine sd_verlet_ciccotti_setup
3632 !-----------------------------------------------------------------------------
3633       subroutine sd_verlet1_ciccotti
3634
3635 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3636 !      implicit real*8 (a-h,o-z)
3637       use energy_data
3638 !      include 'DIMENSIONS'
3639 #ifdef MPI
3640       include 'mpif.h'
3641 #endif
3642 !      include 'COMMON.CONTROL'
3643 !      include 'COMMON.VAR'
3644 !      include 'COMMON.MD'
3645 !#ifndef LANG0
3646 !      include 'COMMON.LANGEVIN'
3647 !#else
3648 !      include 'COMMON.LANGEVIN.lang0'
3649 !#endif
3650 !      include 'COMMON.CHAIN'
3651 !      include 'COMMON.DERIV'
3652 !      include 'COMMON.GEO'
3653 !      include 'COMMON.LOCAL'
3654 !      include 'COMMON.INTERACT'
3655 !      include 'COMMON.IOUNITS'
3656 !      include 'COMMON.NAMES'
3657 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3658 !el      common /stochcalc/ stochforcvec
3659       logical :: lprn = .false.
3660       real(kind=8) :: ddt1,ddt2
3661       integer :: i,j,ind,inres
3662 !      write (iout,*) "dc_old"
3663 !      do i=0,nres
3664 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3665 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3666 !      enddo
3667       do j=1,3
3668         dc_work(j)=dc_old(j,0)
3669         d_t_work(j)=d_t_old(j,0)
3670         d_a_work(j)=d_a_old(j,0)
3671       enddo
3672       ind=3
3673       do i=nnt,nct-1
3674         do j=1,3
3675           dc_work(ind+j)=dc_old(j,i)
3676           d_t_work(ind+j)=d_t_old(j,i)
3677           d_a_work(ind+j)=d_a_old(j,i)
3678         enddo
3679         ind=ind+3
3680       enddo
3681       do i=nnt,nct
3682         if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3683           do j=1,3
3684             dc_work(ind+j)=dc_old(j,i+nres)
3685             d_t_work(ind+j)=d_t_old(j,i+nres)
3686             d_a_work(ind+j)=d_a_old(j,i+nres)
3687           enddo
3688           ind=ind+3
3689         endif
3690       enddo
3691
3692 #ifndef LANG0
3693       if (lprn) then
3694       write (iout,*) &
3695         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3696         " vrand_mat2"
3697       do i=1,dimen
3698         do j=1,dimen
3699           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3700             vfric_mat(i,j),afric_mat(i,j),&
3701             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3702         enddo
3703       enddo
3704       endif
3705       do i=1,dimen
3706         ddt1=0.0d0
3707         ddt2=0.0d0
3708         do j=1,dimen
3709           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3710             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3711           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3712           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3713         enddo
3714         d_t_work_new(i)=ddt1+0.5d0*ddt2
3715         d_t_work(i)=ddt1+ddt2
3716       enddo
3717 #endif
3718       do j=1,3
3719         dc(j,0)=dc_work(j)
3720         d_t(j,0)=d_t_work(j)
3721       enddo
3722       ind=3
3723       do i=nnt,nct-1
3724         do j=1,3
3725           dc(j,i)=dc_work(ind+j)
3726           d_t(j,i)=d_t_work(ind+j)
3727         enddo
3728         ind=ind+3
3729       enddo
3730       do i=nnt,nct
3731          mnum=molnum(i)
3732 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3733          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3734           .and.(mnum.ne.5)) then
3735           inres=i+nres
3736           do j=1,3
3737             dc(j,inres)=dc_work(ind+j)
3738             d_t(j,inres)=d_t_work(ind+j)
3739           enddo
3740           ind=ind+3
3741         endif
3742       enddo
3743       return
3744       end subroutine sd_verlet1_ciccotti
3745 !-----------------------------------------------------------------------------
3746       subroutine sd_verlet2_ciccotti
3747
3748 !  Calculating the adjusted velocities for accelerations
3749       use energy_data
3750 !      implicit real*8 (a-h,o-z)
3751 !      include 'DIMENSIONS'
3752 !      include 'COMMON.CONTROL'
3753 !      include 'COMMON.VAR'
3754 !      include 'COMMON.MD'
3755 !#ifndef LANG0
3756 !      include 'COMMON.LANGEVIN'
3757 !#else
3758 !      include 'COMMON.LANGEVIN.lang0'
3759 !#endif
3760 !      include 'COMMON.CHAIN'
3761 !      include 'COMMON.DERIV'
3762 !      include 'COMMON.GEO'
3763 !      include 'COMMON.LOCAL'
3764 !      include 'COMMON.INTERACT'
3765 !      include 'COMMON.IOUNITS'
3766 !      include 'COMMON.NAMES'
3767 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3768        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3769 !el      common /stochcalc/ stochforcvec
3770       real(kind=8) :: ddt1,ddt2
3771       integer :: i,j,ind,inres
3772 !
3773 ! Compute the stochastic forces which contribute to velocity change
3774 !
3775       call stochastic_force(stochforcvecV)
3776 #ifndef LANG0
3777       do i=1,dimen
3778         ddt1=0.0d0
3779         ddt2=0.0d0
3780         do j=1,dimen
3781
3782           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3783 !          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3784           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3785         enddo
3786         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3787       enddo
3788 #endif
3789       do j=1,3
3790         d_t(j,0)=d_t_work(j)
3791       enddo
3792       ind=3
3793       do i=nnt,nct-1
3794         do j=1,3
3795           d_t(j,i)=d_t_work(ind+j)
3796         enddo
3797         ind=ind+3
3798       enddo
3799       do i=nnt,nct
3800          mnum=molnum(i)
3801          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3802           .and.(mnum.ne.5))
3803 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3804           inres=i+nres
3805           do j=1,3
3806             d_t(j,inres)=d_t_work(ind+j)
3807           enddo
3808           ind=ind+3
3809         endif
3810       enddo
3811       return
3812       end subroutine sd_verlet2_ciccotti
3813 #endif
3814 !-----------------------------------------------------------------------------
3815 ! moments.f
3816 !-----------------------------------------------------------------------------
3817       subroutine inertia_tensor
3818
3819 ! Calculating the intertia tensor for the entire protein in order to
3820 ! remove the perpendicular components of velocity matrix which cause
3821 ! the molecule to rotate.
3822       use comm_gucio
3823       use energy_data
3824 !       implicit real*8 (a-h,o-z)
3825 !       include 'DIMENSIONS'
3826 !       include 'COMMON.CONTROL'
3827 !       include 'COMMON.VAR'
3828 !       include 'COMMON.MD'
3829 !       include 'COMMON.CHAIN'
3830 !       include 'COMMON.DERIV'
3831 !       include 'COMMON.GEO'
3832 !       include 'COMMON.LOCAL'
3833 !       include 'COMMON.INTERACT'
3834 !       include 'COMMON.IOUNITS'
3835 !       include 'COMMON.NAMES'
3836
3837       real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3838       real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3839       real(kind=8) :: M_SC,mag,mag2,M_PEP
3840       real(kind=8),dimension(3,0:nres) :: vpp   !(3,0:MAXRES)
3841       real(kind=8),dimension(3) :: vs_p,pp,incr,v
3842       real(kind=8),dimension(3,3) :: pr1,pr2
3843
3844 !el      common /gucio/ cm
3845       integer :: iti,inres,i,j,k,mnum
3846         do i=1,3
3847            do j=1,3
3848              Im(i,j)=0.0d0
3849              pr1(i,j)=0.0d0
3850              pr2(i,j)=0.0d0
3851            enddo
3852            L(i)=0.0d0
3853            cm(i)=0.0d0
3854            vrot(i)=0.0d0
3855         enddo
3856 !   calculating the center of the mass of the protein
3857         M_PEP=0.0d0
3858         do i=nnt,nct-1
3859           mnum=molnum(i)
3860           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3861           if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3862           M_PEP=M_PEP+mp(mnum)
3863           do j=1,3
3864             cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3865           enddo
3866         enddo
3867 !        do j=1,3
3868 !         cm(j)=mp(1)*cm(j)
3869 !        enddo
3870         M_SC=0.0d0
3871         do i=nnt,nct
3872            mnum=molnum(i)
3873            if (mnum.eq.5) cycle
3874            iti=iabs(itype(i,mnum))
3875            M_SC=M_SC+msc(iabs(iti),mnum)
3876            inres=i+nres
3877            do j=1,3
3878             cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3879            enddo
3880         enddo
3881         do j=1,3
3882           cm(j)=cm(j)/(M_SC+M_PEP)
3883         enddo
3884
3885         do i=nnt,nct-1
3886            mnum=molnum(i)
3887           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3888           do j=1,3
3889             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3890           enddo
3891           Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3892           Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3893           Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3894           Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3895           Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3896           Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3897         enddo
3898
3899         do i=nnt,nct
3900            mnum=molnum(i)
3901            iti=iabs(itype(i,mnum))
3902            if (mnum.eq.5) cycle
3903            inres=i+nres
3904            do j=1,3
3905              pr(j)=c(j,inres)-cm(j)
3906            enddo
3907           Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3908           Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3909           Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3910           Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3911           Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3912           Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3913         enddo
3914
3915         do i=nnt,nct-1
3916            mnum=molnum(i)
3917           Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3918           vbld(i+1)*vbld(i+1)*0.25d0
3919           Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3920           vbld(i+1)*vbld(i+1)*0.25d0
3921           Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3922           vbld(i+1)*vbld(i+1)*0.25d0
3923           Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3924           vbld(i+1)*vbld(i+1)*0.25d0
3925           Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3926           vbld(i+1)*vbld(i+1)*0.25d0
3927           Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3928           vbld(i+1)*vbld(i+1)*0.25d0
3929         enddo
3930
3931
3932         do i=nnt,nct
3933               mnum=molnum(i)
3934 !         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3935          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3936           .and.(mnum.ne.5)) then
3937            iti=iabs(itype(i,mnum))
3938            inres=i+nres
3939           Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3940           dc_norm(1,inres))*vbld(inres)*vbld(inres)
3941           Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3942           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3943           Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3944           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3945           Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3946           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3947           Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3948           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3949           Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3950           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3951          endif
3952         enddo
3953
3954         call angmom(cm,L)
3955 !        write(iout,*) "The angular momentum before adjustment:"
3956 !        write(iout,*) (L(j),j=1,3)
3957
3958         Im(2,1)=Im(1,2)
3959         Im(3,1)=Im(1,3)
3960         Im(3,2)=Im(2,3)
3961
3962 !  Copying the Im matrix for the djacob subroutine
3963         do i=1,3
3964           do j=1,3
3965             Imcp(i,j)=Im(i,j)
3966             Id(i,j)=0.0d0
3967           enddo
3968         enddo
3969
3970 !   Finding the eigenvectors and eignvalues of the inertia tensor
3971        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3972 !       write (iout,*) "Eigenvalues & Eigenvectors"
3973 !       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3974 !       write (iout,*)
3975 !       do i=1,3
3976 !         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3977 !       enddo
3978 !   Constructing the diagonalized matrix
3979        do i=1,3
3980          if (dabs(eigval(i)).gt.1.0d-15) then
3981            Id(i,i)=1.0d0/eigval(i)
3982          else
3983            Id(i,i)=0.0d0
3984          endif
3985        enddo
3986         do i=1,3
3987            do j=1,3
3988               Imcp(i,j)=eigvec(j,i)
3989            enddo
3990         enddo
3991         do i=1,3
3992            do j=1,3
3993               do k=1,3
3994                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3995               enddo
3996            enddo
3997         enddo
3998         do i=1,3
3999            do j=1,3
4000               do k=1,3
4001                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4002               enddo
4003            enddo
4004         enddo
4005 !  Calculating the total rotational velocity of the molecule
4006        do i=1,3
4007          do j=1,3
4008            vrot(i)=vrot(i)+pr2(i,j)*L(j)
4009          enddo
4010        enddo
4011 !   Resetting the velocities
4012        do i=nnt,nct-1
4013          call vecpr(vrot(1),dc(1,i),vp)
4014          do j=1,3
4015            d_t(j,i)=d_t(j,i)-vp(j)
4016           enddo
4017         enddo
4018         do i=nnt,nct
4019               mnum=molnum(i)
4020          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4021           .and.(mnum.ne.5)) then
4022 !         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4023 !       if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4024            inres=i+nres
4025            call vecpr(vrot(1),dc(1,inres),vp)
4026            do j=1,3
4027              d_t(j,inres)=d_t(j,inres)-vp(j)
4028            enddo
4029         endif
4030        enddo
4031        call angmom(cm,L)
4032 !       write(iout,*) "The angular momentum after adjustment:"
4033 !       write(iout,*) (L(j),j=1,3)
4034
4035       return
4036       end subroutine inertia_tensor
4037 !-----------------------------------------------------------------------------
4038       subroutine angmom(cm,L)
4039
4040       use energy_data
4041 !       implicit real*8 (a-h,o-z)
4042 !       include 'DIMENSIONS'
4043 !       include 'COMMON.CONTROL'
4044 !       include 'COMMON.VAR'
4045 !       include 'COMMON.MD'
4046 !       include 'COMMON.CHAIN'
4047 !       include 'COMMON.DERIV'
4048 !       include 'COMMON.GEO'
4049 !       include 'COMMON.LOCAL'
4050 !       include 'COMMON.INTERACT'
4051 !       include 'COMMON.IOUNITS'
4052 !       include 'COMMON.NAMES'
4053       real(kind=8) :: mscab
4054       real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4055       integer :: iti,inres,i,j,mnum
4056 !  Calculate the angular momentum
4057        do j=1,3
4058           L(j)=0.0d0
4059        enddo
4060        do j=1,3
4061           incr(j)=d_t(j,0)
4062        enddo
4063        do i=nnt,nct-1
4064           mnum=molnum(i)
4065           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4066           do j=1,3
4067             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4068           enddo
4069           do j=1,3
4070             v(j)=incr(j)+0.5d0*d_t(j,i)
4071           enddo
4072           do j=1,3
4073             incr(j)=incr(j)+d_t(j,i)
4074           enddo
4075           call vecpr(pr(1),v(1),vp)
4076           do j=1,3
4077             L(j)=L(j)+mp(mnum)*vp(j)
4078           enddo
4079           do j=1,3
4080              pr(j)=0.5d0*dc(j,i)
4081              pp(j)=0.5d0*d_t(j,i)
4082           enddo
4083          call vecpr(pr(1),pp(1),vp)
4084          do j=1,3
4085              L(j)=L(j)+Ip(mnum)*vp(j)
4086           enddo
4087         enddo
4088         do j=1,3
4089           incr(j)=d_t(j,0)
4090         enddo
4091         do i=nnt,nct
4092           mnum=molnum(i)
4093          iti=iabs(itype(i,mnum))
4094          inres=i+nres
4095         if (mnum.eq.5) then
4096          mscab=0.0d0
4097         else
4098          mscab=msc(iti,mnum)
4099         endif
4100          do j=1,3
4101            pr(j)=c(j,inres)-cm(j)
4102          enddo
4103          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4104           .and.(mnum.ne.5)) then
4105            do j=1,3
4106              v(j)=incr(j)+d_t(j,inres)
4107            enddo
4108          else
4109            do j=1,3
4110              v(j)=incr(j)
4111            enddo
4112          endif
4113          call vecpr(pr(1),v(1),vp)
4114 !         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4115 !           " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4116          do j=1,3
4117             L(j)=L(j)+mscab*vp(j)
4118          enddo
4119 !         write (iout,*) "L",(l(j),j=1,3)
4120          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4121           .and.(mnum.ne.5)) then
4122            do j=1,3
4123             v(j)=incr(j)+d_t(j,inres)
4124            enddo
4125            call vecpr(dc(1,inres),d_t(1,inres),vp)
4126            do j=1,3
4127              L(j)=L(j)+Isc(iti,mnum)*vp(j)
4128           enddo
4129          endif
4130          do j=1,3
4131              incr(j)=incr(j)+d_t(j,i)
4132          enddo
4133        enddo
4134       return
4135       end subroutine angmom
4136 !-----------------------------------------------------------------------------
4137       subroutine vcm_vel(vcm)
4138
4139       use energy_data
4140 !       implicit real*8 (a-h,o-z)
4141 !       include 'DIMENSIONS'
4142 !       include 'COMMON.VAR'
4143 !       include 'COMMON.MD'
4144 !       include 'COMMON.CHAIN'
4145 !       include 'COMMON.DERIV'
4146 !       include 'COMMON.GEO'
4147 !       include 'COMMON.LOCAL'
4148 !       include 'COMMON.INTERACT'
4149 !       include 'COMMON.IOUNITS'
4150        real(kind=8),dimension(3) :: vcm,vv
4151        real(kind=8) :: summas,amas
4152        integer :: i,j,mnum
4153
4154        do j=1,3
4155          vcm(j)=0.0d0
4156          vv(j)=d_t(j,0)
4157        enddo
4158        summas=0.0d0
4159        do i=nnt,nct
4160          mnum=molnum(i)
4161          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4162          if (i.lt.nct) then
4163            summas=summas+mp(mnum)
4164            do j=1,3
4165              vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4166            enddo
4167          endif
4168          if (mnum.ne.5) then
4169          amas=msc(iabs(itype(i,mnum)),mnum)
4170          else
4171          amas=0.0d0
4172          endif
4173          summas=summas+amas
4174          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4175           .and.(mnum.ne.5)) then
4176            do j=1,3
4177              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4178            enddo
4179          else
4180            do j=1,3
4181              vcm(j)=vcm(j)+amas*vv(j)
4182            enddo
4183          endif
4184          do j=1,3
4185            vv(j)=vv(j)+d_t(j,i)
4186          enddo
4187        enddo
4188 !       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4189        do j=1,3
4190          vcm(j)=vcm(j)/summas
4191        enddo
4192       return
4193       end subroutine vcm_vel
4194 !-----------------------------------------------------------------------------
4195 ! rattle.F
4196 !-----------------------------------------------------------------------------
4197       subroutine rattle1
4198 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4199 ! AL 9/24/04
4200       use comm_przech
4201       use energy_data
4202 !      implicit real*8 (a-h,o-z)
4203 !      include 'DIMENSIONS'
4204 #ifdef RATTLE
4205 !      include 'COMMON.CONTROL'
4206 !      include 'COMMON.VAR'
4207 !      include 'COMMON.MD'
4208 !#ifndef LANG0
4209 !      include 'COMMON.LANGEVIN'
4210 !#else
4211 !      include 'COMMON.LANGEVIN.lang0'
4212 !#endif
4213 !      include 'COMMON.CHAIN'
4214 !      include 'COMMON.DERIV'
4215 !      include 'COMMON.GEO'
4216 !      include 'COMMON.LOCAL'
4217 !      include 'COMMON.INTERACT'
4218 !      include 'COMMON.IOUNITS'
4219 !      include 'COMMON.NAMES'
4220 !      include 'COMMON.TIME1'
4221 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4222 !el       gdc(3,2*nres,2*nres)
4223       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4224 !el      real(kind=8) :: Cmat(2*nres,2*nres)
4225       real(kind=8) :: x(2*nres),xcorr(3,2*nres)         !maxres2=2*maxres
4226 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4227 !el      common /przechowalnia/ nbond
4228       integer :: max_rattle = 5
4229       logical :: lprn = .false., lprn1 = .false., not_done
4230       real(kind=8) :: tol_rattle = 1.0d-5
4231
4232       integer :: ii,i,j,jj,l,ind,ind1,nres2
4233       nres2=2*nres
4234
4235 !el /common/ przechowalnia
4236
4237       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4238       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4239       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4240 !el--------
4241       if (lprn) write (iout,*) "RATTLE1"
4242       nbond=nct-nnt
4243       do i=nnt,nct
4244        mnum=molnum(i)
4245          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4246           .and.(mnum.ne.5)) nbond=nbond+1
4247       enddo
4248 ! Make a folded form of the Ginv-matrix
4249       ind=0
4250       ii=0
4251       do i=nnt,nct-1
4252         ii=ii+1
4253         do j=1,3
4254           ind=ind+1
4255           ind1=0
4256           jj=0
4257           do k=nnt,nct-1
4258             jj=jj+1
4259             do l=1,3
4260               ind1=ind1+1
4261               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4262             enddo
4263           enddo
4264           do k=nnt,nct
4265           mnum=molnum(k)
4266          if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4267           .and.(mnum.ne.5)) then
4268               jj=jj+1
4269               do l=1,3
4270                 ind1=ind1+1
4271                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4272               enddo
4273             endif
4274           enddo
4275         enddo
4276       enddo
4277       do i=nnt,nct
4278          mnum=molnum(i)
4279          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4280           .and.(mnum.ne.5))
4281           ii=ii+1
4282           do j=1,3
4283             ind=ind+1
4284             ind1=0
4285             jj=0
4286             do k=nnt,nct-1
4287               jj=jj+1
4288               do l=1,3
4289                 ind1=ind1+1
4290                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4291               enddo
4292             enddo
4293             do k=nnt,nct
4294               if (itype(k,1).ne.10) then
4295                 jj=jj+1
4296                 do l=1,3
4297                   ind1=ind1+1
4298                   if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4299                 enddo
4300               endif
4301             enddo
4302           enddo
4303         endif
4304       enddo
4305       if (lprn1) then
4306         write (iout,*) "Matrix GGinv"
4307         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4308       endif
4309       not_done=.true.
4310       iter=0
4311       do while (not_done)
4312       iter=iter+1
4313       if (iter.gt.max_rattle) then
4314         write (iout,*) "Error - too many iterations in RATTLE."
4315         stop
4316       endif
4317 ! Calculate the matrix C = GG**(-1) dC_old o dC
4318       ind1=0
4319       do i=nnt,nct-1
4320         ind1=ind1+1
4321         do j=1,3
4322           dC_uncor(j,ind1)=dC(j,i)
4323         enddo
4324       enddo
4325       do i=nnt,nct
4326         if (itype(i,1).ne.10) then
4327           ind1=ind1+1
4328           do j=1,3
4329             dC_uncor(j,ind1)=dC(j,i+nres)
4330           enddo
4331         endif
4332       enddo
4333       do i=1,nbond
4334         ind=0
4335         do k=nnt,nct-1
4336           ind=ind+1
4337           do j=1,3
4338             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4339           enddo
4340         enddo
4341         do k=nnt,nct
4342           if (itype(k,1).ne.10) then
4343             ind=ind+1
4344             do j=1,3
4345               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4346             enddo
4347           endif
4348         enddo
4349       enddo
4350 ! Calculate deviations from standard virtual-bond lengths
4351       ind=0
4352       do i=nnt,nct-1
4353         ind=ind+1
4354         x(ind)=vbld(i+1)**2-vbl**2
4355       enddo
4356       do i=nnt,nct
4357         if (itype(i,1).ne.10) then
4358           ind=ind+1
4359           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4360         endif
4361       enddo
4362       if (lprn) then
4363         write (iout,*) "Coordinates and violations"
4364         do i=1,nbond
4365           write(iout,'(i5,3f10.5,5x,e15.5)') &
4366            i,(dC_uncor(j,i),j=1,3),x(i)
4367         enddo
4368         write (iout,*) "Velocities and violations"
4369         ind=0
4370         do i=nnt,nct-1
4371           ind=ind+1
4372           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4373            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4374         enddo
4375         do i=nnt,nct
4376           mnum=molnum(i)
4377          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4378           .and.(mnum.ne.5)) then
4379
4380             ind=ind+1
4381             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4382              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4383              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4384           endif
4385         enddo
4386 !        write (iout,*) "gdc"
4387 !        do i=1,nbond
4388 !          write (iout,*) "i",i
4389 !          do j=1,nbond
4390 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4391 !          enddo
4392 !        enddo
4393       endif
4394       xmax=dabs(x(1))
4395       do i=2,nbond
4396         if (dabs(x(i)).gt.xmax) then
4397           xmax=dabs(x(i))
4398         endif
4399       enddo
4400       if (xmax.lt.tol_rattle) then
4401         not_done=.false.
4402         goto 100
4403       endif
4404 ! Calculate the matrix of the system of equations
4405       do i=1,nbond
4406         do j=1,nbond
4407           Cmat(i,j)=0.0d0
4408           do k=1,3
4409             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4410           enddo
4411         enddo
4412       enddo
4413       if (lprn1) then
4414         write (iout,*) "Matrix Cmat"
4415         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4416       endif
4417       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4418 ! Add constraint term to positions
4419       ind=0
4420       do i=nnt,nct-1
4421         ind=ind+1
4422         do j=1,3
4423           xx=0.0d0
4424           do ii=1,nbond
4425             xx = xx+x(ii)*gdc(j,ind,ii)
4426           enddo
4427           xx=0.5d0*xx
4428           dC(j,i)=dC(j,i)-xx
4429           d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4430         enddo
4431       enddo
4432       do i=nnt,nct
4433         if (itype(i,1).ne.10) then
4434           ind=ind+1
4435           do j=1,3
4436             xx=0.0d0
4437             do ii=1,nbond
4438               xx = xx+x(ii)*gdc(j,ind,ii)
4439             enddo
4440             xx=0.5d0*xx
4441             dC(j,i+nres)=dC(j,i+nres)-xx
4442             d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4443           enddo
4444         endif
4445       enddo
4446 ! Rebuild the chain using the new coordinates
4447       call chainbuild_cart
4448       if (lprn) then
4449         write (iout,*) "New coordinates, Lagrange multipliers,",&
4450         " and differences between actual and standard bond lengths"
4451         ind=0
4452         do i=nnt,nct-1
4453           ind=ind+1
4454           xx=vbld(i+1)**2-vbl**2
4455           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4456               i,(dC(j,i),j=1,3),x(ind),xx
4457         enddo
4458         do i=nnt,nct
4459          mnum=molnum(i)
4460          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4461           .and.(mnum.ne.5))
4462             ind=ind+1
4463             xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4464             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4465              i,(dC(j,i+nres),j=1,3),x(ind),xx
4466           endif
4467         enddo
4468         write (iout,*) "Velocities and violations"
4469         ind=0
4470         do i=nnt,nct-1
4471           ind=ind+1
4472           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4473            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4474         enddo
4475         do i=nnt,nct
4476           if (itype(i,1).ne.10) then
4477             ind=ind+1
4478             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4479              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4480              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4481           endif
4482         enddo
4483       endif
4484       enddo
4485   100 continue
4486       return
4487    10 write (iout,*) "Error - singularity in solving the system",&
4488        " of equations for Lagrange multipliers."
4489       stop
4490 #else
4491       write (iout,*) &
4492        "RATTLE inactive; use -DRATTLE switch at compile time."
4493       stop
4494 #endif
4495       end subroutine rattle1
4496 !-----------------------------------------------------------------------------
4497       subroutine rattle2
4498 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4499 ! AL 9/24/04
4500       use comm_przech
4501       use energy_data
4502 !      implicit real*8 (a-h,o-z)
4503 !      include 'DIMENSIONS'
4504 #ifdef RATTLE
4505 !      include 'COMMON.CONTROL'
4506 !      include 'COMMON.VAR'
4507 !      include 'COMMON.MD'
4508 !#ifndef LANG0
4509 !      include 'COMMON.LANGEVIN'
4510 !#else
4511 !      include 'COMMON.LANGEVIN.lang0'
4512 !#endif
4513 !      include 'COMMON.CHAIN'
4514 !      include 'COMMON.DERIV'
4515 !      include 'COMMON.GEO'
4516 !      include 'COMMON.LOCAL'
4517 !      include 'COMMON.INTERACT'
4518 !      include 'COMMON.IOUNITS'
4519 !      include 'COMMON.NAMES'
4520 !      include 'COMMON.TIME1'
4521 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4522 !el       gdc(3,2*nres,2*nres)
4523       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4524 !el       Cmat(2*nres,2*nres)
4525       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4526 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4527 !el      common /przechowalnia/ nbond
4528       integer :: max_rattle = 5
4529       logical :: lprn = .false., lprn1 = .false., not_done
4530       real(kind=8) :: tol_rattle = 1.0d-5
4531       integer :: nres2
4532       nres2=2*nres
4533
4534 !el /common/ przechowalnia
4535       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4536       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4537       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4538 !el--------
4539       if (lprn) write (iout,*) "RATTLE2"
4540       if (lprn) write (iout,*) "Velocity correction"
4541 ! Calculate the matrix G dC
4542       do i=1,nbond
4543         ind=0
4544         do k=nnt,nct-1
4545           ind=ind+1
4546           do j=1,3
4547             gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4548           enddo
4549         enddo
4550         do k=nnt,nct
4551          mnum=molnum(i)
4552          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4553           .and.(mnum.ne.5)) then
4554             ind=ind+1
4555             do j=1,3
4556               gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4557             enddo
4558           endif
4559         enddo
4560       enddo
4561 !      if (lprn) then
4562 !        write (iout,*) "gdc"
4563 !        do i=1,nbond
4564 !          write (iout,*) "i",i
4565 !          do j=1,nbond
4566 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4567 !          enddo
4568 !        enddo
4569 !      endif
4570 ! Calculate the matrix of the system of equations
4571       ind=0
4572       do i=nnt,nct-1
4573         ind=ind+1
4574         do j=1,nbond
4575           Cmat(ind,j)=0.0d0
4576           do k=1,3
4577             Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4578           enddo
4579         enddo
4580       enddo
4581       do i=nnt,nct
4582          mnum=molnum(i)
4583          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4584           .and.(mnum.ne.5)) then
4585           ind=ind+1
4586           do j=1,nbond
4587             Cmat(ind,j)=0.0d0
4588             do k=1,3
4589               Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4590             enddo
4591           enddo
4592         endif
4593       enddo
4594 ! Calculate the scalar product dC o d_t_new
4595       ind=0
4596       do i=nnt,nct-1
4597         ind=ind+1
4598         x(ind)=scalar(d_t(1,i),dC(1,i))
4599       enddo
4600       do i=nnt,nct
4601          mnum=molnum(i)
4602          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4603           .and.(mnum.ne.5)) then
4604           ind=ind+1
4605           x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4606         endif
4607       enddo
4608       if (lprn) then
4609         write (iout,*) "Velocities and violations"
4610         ind=0
4611         do i=nnt,nct-1
4612           ind=ind+1
4613           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4614            i,ind,(d_t(j,i),j=1,3),x(ind)
4615         enddo
4616         do i=nnt,nct
4617          mnum=molnum(i)
4618          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4619           .and.(mnum.ne.5)) then
4620             ind=ind+1
4621             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4622              i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4623           endif
4624         enddo
4625       endif
4626       xmax=dabs(x(1))
4627       do i=2,nbond
4628         if (dabs(x(i)).gt.xmax) then
4629           xmax=dabs(x(i))
4630         endif
4631       enddo
4632       if (xmax.lt.tol_rattle) then
4633         not_done=.false.
4634         goto 100
4635       endif
4636       if (lprn1) then
4637         write (iout,*) "Matrix Cmat"
4638         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4639       endif
4640       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4641 ! Add constraint term to velocities
4642       ind=0
4643       do i=nnt,nct-1
4644         ind=ind+1
4645         do j=1,3
4646           xx=0.0d0
4647           do ii=1,nbond
4648             xx = xx+x(ii)*gdc(j,ind,ii)
4649           enddo
4650           d_t(j,i)=d_t(j,i)-xx
4651         enddo
4652       enddo
4653       do i=nnt,nct
4654          mnum=molnum(i)
4655          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4656           .and.(mnum.ne.5)) then
4657           ind=ind+1
4658           do j=1,3
4659             xx=0.0d0
4660             do ii=1,nbond
4661               xx = xx+x(ii)*gdc(j,ind,ii)
4662             enddo
4663             d_t(j,i+nres)=d_t(j,i+nres)-xx
4664           enddo
4665         endif
4666       enddo
4667       if (lprn) then
4668         write (iout,*) &
4669           "New velocities, Lagrange multipliers violations"
4670         ind=0
4671         do i=nnt,nct-1
4672           ind=ind+1
4673           if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4674              i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4675         enddo
4676         do i=nnt,nct
4677          mnum=molnum(i)
4678          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4679           .and.(mnum.ne.5))
4680             ind=ind+1
4681             write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4682               i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4683               scalar(d_t(1,i+nres),dC(1,i+nres))
4684           endif
4685         enddo
4686       endif
4687   100 continue
4688       return
4689    10 write (iout,*) "Error - singularity in solving the system",&
4690        " of equations for Lagrange multipliers."
4691       stop
4692 #else
4693       write (iout,*) &
4694        "RATTLE inactive; use -DRATTLE option at compile time."
4695       stop
4696 #endif
4697       end subroutine rattle2
4698 !-----------------------------------------------------------------------------
4699       subroutine rattle_brown
4700 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4701 ! AL 9/24/04
4702       use comm_przech
4703       use energy_data
4704 !      implicit real*8 (a-h,o-z)
4705 !      include 'DIMENSIONS'
4706 #ifdef RATTLE
4707 !      include 'COMMON.CONTROL'
4708 !      include 'COMMON.VAR'
4709 !      include 'COMMON.MD'
4710 !#ifndef LANG0
4711 !      include 'COMMON.LANGEVIN'
4712 !#else
4713 !      include 'COMMON.LANGEVIN.lang0'
4714 !#endif
4715 !      include 'COMMON.CHAIN'
4716 !      include 'COMMON.DERIV'
4717 !      include 'COMMON.GEO'
4718 !      include 'COMMON.LOCAL'
4719 !      include 'COMMON.INTERACT'
4720 !      include 'COMMON.IOUNITS'
4721 !      include 'COMMON.NAMES'
4722 !      include 'COMMON.TIME1'
4723 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4724 !el       gdc(3,2*nres,2*nres)
4725       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4726 !el      real(kind=8) :: Cmat(2*nres,2*nres)
4727       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4728 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4729 !el      common /przechowalnia/ nbond
4730       integer :: max_rattle = 5
4731       logical :: lprn = .false., lprn1 = .false., not_done
4732       real(kind=8) :: tol_rattle = 1.0d-5
4733       integer :: nres2
4734       nres2=2*nres
4735
4736 !el /common/ przechowalnia
4737       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4738       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4739       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4740 !el--------
4741
4742       if (lprn) write (iout,*) "RATTLE_BROWN"
4743       nbond=nct-nnt
4744       do i=nnt,nct
4745         if (itype(i,1).ne.10) nbond=nbond+1
4746       enddo
4747 ! Make a folded form of the Ginv-matrix
4748       ind=0
4749       ii=0
4750       do i=nnt,nct-1
4751         ii=ii+1
4752         do j=1,3
4753           ind=ind+1
4754           ind1=0
4755           jj=0
4756           do k=nnt,nct-1
4757             jj=jj+1
4758             do l=1,3
4759               ind1=ind1+1
4760               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4761             enddo
4762           enddo
4763           do k=nnt,nct
4764             if (itype(k,1).ne.10) then
4765               jj=jj+1
4766               do l=1,3
4767                 ind1=ind1+1
4768                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4769               enddo
4770             endif
4771           enddo
4772         enddo
4773       enddo
4774       do i=nnt,nct
4775         if (itype(i,1).ne.10) then
4776           ii=ii+1
4777           do j=1,3
4778             ind=ind+1
4779             ind1=0
4780             jj=0
4781             do k=nnt,nct-1
4782               jj=jj+1
4783               do l=1,3
4784                 ind1=ind1+1
4785                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4786               enddo
4787             enddo
4788             do k=nnt,nct
4789               if (itype(k,1).ne.10) then
4790                 jj=jj+1
4791                 do l=1,3
4792                   ind1=ind1+1
4793                   if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4794                 enddo
4795               endif
4796             enddo
4797           enddo
4798         endif
4799       enddo
4800       if (lprn1) then
4801         write (iout,*) "Matrix GGinv"
4802         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4803       endif
4804       not_done=.true.
4805       iter=0
4806       do while (not_done)
4807       iter=iter+1
4808       if (iter.gt.max_rattle) then
4809         write (iout,*) "Error - too many iterations in RATTLE."
4810         stop
4811       endif
4812 ! Calculate the matrix C = GG**(-1) dC_old o dC
4813       ind1=0
4814       do i=nnt,nct-1
4815         ind1=ind1+1
4816         do j=1,3
4817           dC_uncor(j,ind1)=dC(j,i)
4818         enddo
4819       enddo
4820       do i=nnt,nct
4821         if (itype(i,1).ne.10) then
4822           ind1=ind1+1
4823           do j=1,3
4824             dC_uncor(j,ind1)=dC(j,i+nres)
4825           enddo
4826         endif
4827       enddo
4828       do i=1,nbond
4829         ind=0
4830         do k=nnt,nct-1
4831           ind=ind+1
4832           do j=1,3
4833             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4834           enddo
4835         enddo
4836         do k=nnt,nct
4837           if (itype(k,1).ne.10) then
4838             ind=ind+1
4839             do j=1,3
4840               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4841             enddo
4842           endif
4843         enddo
4844       enddo
4845 ! Calculate deviations from standard virtual-bond lengths
4846       ind=0
4847       do i=nnt,nct-1
4848         ind=ind+1
4849         x(ind)=vbld(i+1)**2-vbl**2
4850       enddo
4851       do i=nnt,nct
4852         if (itype(i,1).ne.10) then
4853           ind=ind+1
4854           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4855         endif
4856       enddo
4857       if (lprn) then
4858         write (iout,*) "Coordinates and violations"
4859         do i=1,nbond
4860           write(iout,'(i5,3f10.5,5x,e15.5)') &
4861            i,(dC_uncor(j,i),j=1,3),x(i)
4862         enddo
4863         write (iout,*) "Velocities and violations"
4864         ind=0
4865         do i=nnt,nct-1
4866           ind=ind+1
4867           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4868            i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4869         enddo
4870         do i=nnt,nct
4871           if (itype(i,1).ne.10) then
4872             ind=ind+1
4873             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4874              i+nres,ind,(d_t(j,i+nres),j=1,3),&
4875              scalar(d_t(1,i+nres),dC_old(1,i+nres))
4876           endif
4877         enddo
4878         write (iout,*) "gdc"
4879         do i=1,nbond
4880           write (iout,*) "i",i
4881           do j=1,nbond
4882             write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4883           enddo
4884         enddo
4885       endif
4886       xmax=dabs(x(1))
4887       do i=2,nbond
4888         if (dabs(x(i)).gt.xmax) then
4889           xmax=dabs(x(i))
4890         endif
4891       enddo
4892       if (xmax.lt.tol_rattle) then
4893         not_done=.false.
4894         goto 100
4895       endif
4896 ! Calculate the matrix of the system of equations
4897       do i=1,nbond
4898         do j=1,nbond
4899           Cmat(i,j)=0.0d0
4900           do k=1,3
4901             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4902           enddo
4903         enddo
4904       enddo
4905       if (lprn1) then
4906         write (iout,*) "Matrix Cmat"
4907         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4908       endif
4909       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4910 ! Add constraint term to positions
4911       ind=0
4912       do i=nnt,nct-1
4913         ind=ind+1
4914         do j=1,3
4915           xx=0.0d0
4916           do ii=1,nbond
4917             xx = xx+x(ii)*gdc(j,ind,ii)
4918           enddo
4919           xx=-0.5d0*xx
4920           d_t(j,i)=d_t(j,i)+xx/d_time
4921           dC(j,i)=dC(j,i)+xx
4922         enddo
4923       enddo
4924       do i=nnt,nct
4925         if (itype(i,1).ne.10) then
4926           ind=ind+1
4927           do j=1,3
4928             xx=0.0d0
4929             do ii=1,nbond
4930               xx = xx+x(ii)*gdc(j,ind,ii)
4931             enddo
4932             xx=-0.5d0*xx
4933             d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4934             dC(j,i+nres)=dC(j,i+nres)+xx
4935           enddo
4936         endif
4937       enddo
4938 ! Rebuild the chain using the new coordinates
4939       call chainbuild_cart
4940       if (lprn) then
4941         write (iout,*) "New coordinates, Lagrange multipliers,",&
4942         " and differences between actual and standard bond lengths"
4943         ind=0
4944         do i=nnt,nct-1
4945           ind=ind+1
4946           xx=vbld(i+1)**2-vbl**2
4947           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4948               i,(dC(j,i),j=1,3),x(ind),xx
4949         enddo
4950         do i=nnt,nct
4951           if (itype(i,1).ne.10) then
4952             ind=ind+1
4953             xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4954             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4955              i,(dC(j,i+nres),j=1,3),x(ind),xx
4956           endif
4957         enddo
4958         write (iout,*) "Velocities and violations"
4959         ind=0
4960         do i=nnt,nct-1
4961           ind=ind+1
4962           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4963            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4964         enddo
4965         do i=nnt,nct
4966           if (itype(i,1).ne.10) then
4967             ind=ind+1
4968             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4969              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4970              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4971           endif
4972         enddo
4973       endif
4974       enddo
4975   100 continue
4976       return
4977    10 write (iout,*) "Error - singularity in solving the system",&
4978        " of equations for Lagrange multipliers."
4979       stop
4980 #else
4981       write (iout,*) &
4982        "RATTLE inactive; use -DRATTLE option at compile time"
4983       stop
4984 #endif
4985       end subroutine rattle_brown
4986 !-----------------------------------------------------------------------------
4987 ! stochfric.F
4988 !-----------------------------------------------------------------------------
4989       subroutine friction_force
4990
4991       use energy_data
4992       use REMD_data
4993       use comm_syfek
4994 !      implicit real*8 (a-h,o-z)
4995 !      include 'DIMENSIONS'
4996 !      include 'COMMON.VAR'
4997 !      include 'COMMON.CHAIN'
4998 !      include 'COMMON.DERIV'
4999 !      include 'COMMON.GEO'
5000 !      include 'COMMON.LOCAL'
5001 !      include 'COMMON.INTERACT'
5002 !      include 'COMMON.MD'
5003 !#ifndef LANG0
5004 !      include 'COMMON.LANGEVIN'
5005 !#else
5006 !      include 'COMMON.LANGEVIN.lang0'
5007 !#endif
5008 !      include 'COMMON.IOUNITS'
5009 !el      real(kind=8),dimension(6*nres) :: gamvec       !(MAXRES6) maxres6=6*maxres
5010 !el      common /syfek/ gamvec
5011       real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5012 !el       ginvfric(2*nres,2*nres)       !maxres2=2*maxres
5013 !el      common /przechowalnia/ ginvfric
5014
5015       logical :: lprn = .false., checkmode = .false.
5016       integer :: i,j,ind,k,nres2,nres6,mnum
5017       nres2=2*nres
5018       nres6=6*nres
5019
5020       if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5021       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5022       do i=0,nres2
5023         do j=1,3
5024           friction(j,i)=0.0d0
5025         enddo
5026       enddo
5027
5028       do j=1,3
5029         d_t_work(j)=d_t(j,0)
5030       enddo
5031       ind=3
5032       do i=nnt,nct-1
5033         do j=1,3
5034           d_t_work(ind+j)=d_t(j,i)
5035         enddo
5036         ind=ind+3
5037       enddo
5038       do i=nnt,nct
5039         mnum=molnum(i)
5040         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5041         .and.(mnum.ne.5)) then
5042           do j=1,3
5043             d_t_work(ind+j)=d_t(j,i+nres)
5044           enddo
5045           ind=ind+3
5046         endif
5047       enddo
5048
5049       call fricmat_mult(d_t_work,fric_work)
5050
5051       if (.not.checkmode) return
5052
5053       if (lprn) then
5054         write (iout,*) "d_t_work and fric_work"
5055         do i=1,3*dimen
5056           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5057         enddo
5058       endif
5059       do j=1,3
5060         friction(j,0)=fric_work(j)
5061       enddo
5062       ind=3
5063       do i=nnt,nct-1
5064         do j=1,3
5065           friction(j,i)=fric_work(ind+j)
5066         enddo
5067         ind=ind+3
5068       enddo
5069       do i=nnt,nct
5070         mnum=molnum(i)
5071         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5072         .and.(mnum.ne.5)) then
5073 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5074           do j=1,3
5075             friction(j,i+nres)=fric_work(ind+j)
5076           enddo
5077           ind=ind+3
5078         endif
5079       enddo
5080       if (lprn) then
5081         write(iout,*) "Friction backbone"
5082         do i=0,nct-1
5083           write(iout,'(i5,3e15.5,5x,3e15.5)') &
5084            i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5085         enddo
5086         write(iout,*) "Friction side chain"
5087         do i=nnt,nct
5088           write(iout,'(i5,3e15.5,5x,3e15.5)') &
5089            i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5090         enddo
5091       endif
5092       if (lprn) then
5093         do j=1,3
5094           vv(j)=d_t(j,0)
5095         enddo
5096         do i=nnt,nct
5097           do j=1,3
5098             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5099             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5100             vv(j)=vv(j)+d_t(j,i)
5101           enddo
5102         enddo
5103         write (iout,*) "vvtot backbone and sidechain"
5104         do i=nnt,nct
5105           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5106            (vvtot(j,i+nres),j=1,3)
5107         enddo
5108         ind=0
5109         do i=nnt,nct-1
5110           do j=1,3
5111             v_work(ind+j)=vvtot(j,i)
5112           enddo
5113           ind=ind+3
5114         enddo
5115         do i=nnt,nct
5116           do j=1,3
5117             v_work(ind+j)=vvtot(j,i+nres)
5118           enddo
5119           ind=ind+3
5120         enddo
5121         write (iout,*) "v_work gamvec and site-based friction forces"
5122         do i=1,dimen1
5123           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5124             gamvec(i)*v_work(i)
5125         enddo
5126 !        do i=1,dimen
5127 !          fric_work1(i)=0.0d0
5128 !          do j=1,dimen1
5129 !            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5130 !          enddo
5131 !        enddo
5132 !        write (iout,*) "fric_work and fric_work1"
5133 !        do i=1,dimen
5134 !          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5135 !        enddo
5136         do i=1,dimen
5137           do j=1,dimen
5138             ginvfric(i,j)=0.0d0
5139             do k=1,dimen
5140               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5141             enddo
5142           enddo
5143         enddo
5144         write (iout,*) "ginvfric"
5145         do i=1,dimen
5146           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5147         enddo
5148         write (iout,*) "symmetry check"
5149         do i=1,dimen
5150           do j=1,i-1
5151             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5152           enddo
5153         enddo
5154       endif
5155       return
5156       end subroutine friction_force
5157 !-----------------------------------------------------------------------------
5158       subroutine setup_fricmat
5159
5160 !     use MPI
5161       use energy_data
5162       use control_data, only:time_Bcast
5163       use control, only:tcpu
5164       use comm_syfek
5165 !      implicit real*8 (a-h,o-z)
5166 #ifdef MPI
5167       use MPI_data
5168       include 'mpif.h'
5169       real(kind=8) :: time00
5170 #endif
5171 !      include 'DIMENSIONS'
5172 !      include 'COMMON.VAR'
5173 !      include 'COMMON.CHAIN'
5174 !      include 'COMMON.DERIV'
5175 !      include 'COMMON.GEO'
5176 !      include 'COMMON.LOCAL'
5177 !      include 'COMMON.INTERACT'
5178 !      include 'COMMON.MD'
5179 !      include 'COMMON.SETUP'
5180 !      include 'COMMON.TIME1'
5181 !      integer licznik /0/
5182 !      save licznik
5183 !#ifndef LANG0
5184 !      include 'COMMON.LANGEVIN'
5185 !#else
5186 !      include 'COMMON.LANGEVIN.lang0'
5187 !#endif
5188 !      include 'COMMON.IOUNITS'
5189       integer :: IERROR
5190       integer :: i,j,ind,ind1,m
5191       logical :: lprn = .false.
5192       real(kind=8) :: dtdi !el ,gamvec(2*nres)
5193 !el      real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5194 !      real(kind=8),allocatable,dimension(:,:) :: fcopy
5195 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5196 !el      common /syfek/ gamvec
5197       real(kind=8) :: work(8*2*nres)
5198       integer :: iwork(2*nres)
5199 !el      common /przechowalnia/ ginvfric,Ghalf,fcopy
5200       integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5201       nres2=2*nres
5202       nres6=6*nres
5203 #ifdef MPI
5204       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5205        if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5206       if (fg_rank.ne.king) goto 10
5207 #endif
5208 !      nres2=2*nres
5209 !      nres6=6*nres
5210
5211       if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5212       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5213        if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5214 !el      allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5215       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5216
5217       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5218 !  Zeroing out fricmat
5219       do i=1,dimen
5220         do j=1,dimen
5221           fricmat(i,j)=0.0d0
5222         enddo
5223       enddo
5224 !  Load the friction coefficients corresponding to peptide groups
5225       ind1=0
5226       do i=nnt,nct-1
5227         mnum=molnum(i)
5228         ind1=ind1+1
5229         gamvec(ind1)=gamp(mnum)
5230       enddo
5231 !  Load the friction coefficients corresponding to side chains
5232       m=nct-nnt
5233       ind=0
5234       do j=1,2
5235       gamsc(ntyp1_molec(j),j)=1.0d0
5236       enddo
5237       do i=nnt,nct
5238         mnum=molnum(i)
5239         ind=ind+1
5240         ii = ind+m
5241         iti=itype(i,mnum)
5242         gamvec(ii)=gamsc(iabs(iti),mnum)
5243       enddo
5244       if (surfarea) call sdarea(gamvec)
5245 !      if (lprn) then
5246 !        write (iout,*) "Matrix A and vector gamma"
5247 !        do i=1,dimen1
5248 !          write (iout,'(i2,$)') i
5249 !          do j=1,dimen
5250 !            write (iout,'(f4.1,$)') A(i,j)
5251 !          enddo
5252 !          write (iout,'(f8.3)') gamvec(i)
5253 !        enddo
5254 !      endif
5255       if (lprn) then
5256         write (iout,*) "Vector gamvec"
5257         do i=1,dimen1
5258           write (iout,'(i5,f10.5)') i, gamvec(i)
5259         enddo
5260       endif
5261
5262 ! The friction matrix
5263       do k=1,dimen
5264        do i=1,dimen
5265          dtdi=0.0d0
5266          do j=1,dimen1
5267            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5268          enddo
5269          fricmat(k,i)=dtdi
5270        enddo
5271       enddo
5272
5273       if (lprn) then
5274         write (iout,'(//a)') "Matrix fricmat"
5275         call matout2(dimen,dimen,nres2,nres2,fricmat)
5276       endif
5277       if (lang.eq.2 .or. lang.eq.3) then
5278 ! Mass-scale the friction matrix if non-direct integration will be performed
5279       do i=1,dimen
5280         do j=1,dimen
5281           Ginvfric(i,j)=0.0d0
5282           do k=1,dimen
5283             do l=1,dimen
5284               Ginvfric(i,j)=Ginvfric(i,j)+ &
5285                 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5286             enddo
5287           enddo
5288         enddo
5289       enddo
5290 ! Diagonalize the friction matrix
5291       ind=0
5292       do i=1,dimen
5293         do j=1,i
5294           ind=ind+1
5295           Ghalf(ind)=Ginvfric(i,j)
5296         enddo
5297       enddo
5298       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5299         ierr,iwork)
5300       if (lprn) then
5301         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5302           " mass-scaled friction matrix"
5303         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5304       endif
5305 ! Precompute matrices for tinker stochastic integrator
5306 #ifndef LANG0
5307       do i=1,dimen
5308         do j=1,dimen
5309           mt1(i,j)=0.0d0
5310           mt2(i,j)=0.0d0
5311           do k=1,dimen
5312             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5313             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5314           enddo
5315           mt3(j,i)=mt1(i,j)
5316         enddo
5317       enddo
5318 #endif
5319       else if (lang.eq.4) then
5320 ! Diagonalize the friction matrix
5321       ind=0
5322       do i=1,dimen
5323         do j=1,i
5324           ind=ind+1
5325           Ghalf(ind)=fricmat(i,j)
5326         enddo
5327       enddo
5328       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5329         ierr,iwork)
5330       if (lprn) then
5331         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5332           " friction matrix"
5333         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5334       endif
5335 ! Determine the number of zero eigenvalues of the friction matrix
5336       nzero=max0(dimen-dimen1,0)
5337 !      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5338 !        nzero=nzero+1
5339 !      enddo
5340       write (iout,*) "Number of zero eigenvalues:",nzero
5341       do i=1,dimen
5342         do j=1,dimen
5343           fricmat(i,j)=0.0d0
5344           do k=nzero+1,dimen
5345             fricmat(i,j)=fricmat(i,j) &
5346               +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5347           enddo
5348         enddo
5349       enddo
5350       if (lprn) then
5351         write (iout,'(//a)') "Generalized inverse of fricmat"
5352         call matout(dimen,dimen,nres6,nres6,fricmat)
5353       endif
5354       endif
5355 #ifdef MPI
5356   10  continue
5357       if (nfgtasks.gt.1) then
5358         if (fg_rank.eq.0) then
5359 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5360 #ifdef MPI
5361           time00=MPI_Wtime()
5362 #else
5363           time00=tcpu()
5364 #endif
5365           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5366 #ifdef MPI
5367           time_Bcast=time_Bcast+MPI_Wtime()-time00
5368 #else
5369           time_Bcast=time_Bcast+tcpu()-time00
5370 #endif
5371 !          print *,"Processor",myrank,
5372 !     &       " BROADCAST iorder in SETUP_FRICMAT"
5373         endif
5374 !       licznik=licznik+1
5375         write (iout,*) "setup_fricmat licznik"!,licznik !sp
5376 #ifdef MPI
5377         time00=MPI_Wtime()
5378 #else
5379         time00=tcpu()
5380 #endif
5381 ! Scatter the friction matrix
5382         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5383           nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5384           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5385 #ifdef TIMING
5386 #ifdef MPI
5387         time_scatter=time_scatter+MPI_Wtime()-time00
5388         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5389 #else
5390         time_scatter=time_scatter+tcpu()-time00
5391         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5392 #endif
5393 #endif
5394         do i=1,dimen
5395           do j=1,2*my_ng_count
5396             fricmat(j,i)=fcopy(i,j)
5397           enddo
5398         enddo
5399 !        write (iout,*) "My chunk of fricmat"
5400 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5401       endif
5402 #endif
5403       return
5404       end subroutine setup_fricmat
5405 !-----------------------------------------------------------------------------
5406       subroutine stochastic_force(stochforcvec)
5407
5408       use energy_data
5409       use random, only:anorm_distr
5410 !      implicit real*8 (a-h,o-z)
5411 !      include 'DIMENSIONS'
5412       use control, only: tcpu
5413       use control_data
5414 #ifdef MPI
5415       include 'mpif.h'
5416 #endif
5417 !      include 'COMMON.VAR'
5418 !      include 'COMMON.CHAIN'
5419 !      include 'COMMON.DERIV'
5420 !      include 'COMMON.GEO'
5421 !      include 'COMMON.LOCAL'
5422 !      include 'COMMON.INTERACT'
5423 !      include 'COMMON.MD'
5424 !      include 'COMMON.TIME1'
5425 !#ifndef LANG0
5426 !      include 'COMMON.LANGEVIN'
5427 !#else
5428 !      include 'COMMON.LANGEVIN.lang0'
5429 !#endif
5430 !      include 'COMMON.IOUNITS'
5431
5432       real(kind=8) :: x,sig,lowb,highb
5433       real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5434       real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5435       real(kind=8) :: time00
5436       logical :: lprn = .false.
5437       integer :: i,j,ind,mnum
5438
5439       do i=0,2*nres
5440         do j=1,3
5441           stochforc(j,i)=0.0d0
5442         enddo
5443       enddo
5444       x=0.0d0
5445
5446 #ifdef MPI
5447       time00=MPI_Wtime()
5448 #else
5449       time00=tcpu()
5450 #endif
5451 ! Compute the stochastic forces acting on bodies. Store in force.
5452       do i=nnt,nct-1
5453         sig=stdforcp(i)
5454         lowb=-5*sig
5455         highb=5*sig
5456         do j=1,3
5457           force(j,i)=anorm_distr(x,sig,lowb,highb)
5458         enddo
5459       enddo
5460       do i=nnt,nct
5461         sig2=stdforcsc(i)
5462         lowb2=-5*sig2
5463         highb2=5*sig2
5464         do j=1,3
5465           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5466         enddo
5467       enddo
5468 #ifdef MPI
5469       time_fsample=time_fsample+MPI_Wtime()-time00
5470 #else
5471       time_fsample=time_fsample+tcpu()-time00
5472 #endif
5473 ! Compute the stochastic forces acting on virtual-bond vectors.
5474       do j=1,3
5475         ff(j)=0.0d0
5476       enddo
5477       do i=nct-1,nnt,-1
5478         do j=1,3
5479           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5480         enddo
5481         do j=1,3
5482           ff(j)=ff(j)+force(j,i)
5483         enddo
5484 !        if (itype(i+1,1).ne.ntyp1) then
5485          mnum=molnum(i)
5486          if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5487           do j=1,3
5488             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5489             ff(j)=ff(j)+force(j,i+nres+1)
5490           enddo
5491         endif
5492       enddo
5493       do j=1,3
5494         stochforc(j,0)=ff(j)+force(j,nnt+nres)
5495       enddo
5496       do i=nnt,nct
5497         mnum=molnum(i)
5498         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5499         .and.(mnum.ne.5)) then
5500 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5501           do j=1,3
5502             stochforc(j,i+nres)=force(j,i+nres)
5503           enddo
5504         endif
5505       enddo
5506
5507       do j=1,3
5508         stochforcvec(j)=stochforc(j,0)
5509       enddo
5510       ind=3
5511       do i=nnt,nct-1
5512         do j=1,3
5513           stochforcvec(ind+j)=stochforc(j,i)
5514         enddo
5515         ind=ind+3
5516       enddo
5517       do i=nnt,nct
5518         mnum=molnum(i)
5519         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5520         .and.(mnum.ne.5)) then
5521 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5522           do j=1,3
5523             stochforcvec(ind+j)=stochforc(j,i+nres)
5524           enddo
5525           ind=ind+3
5526         endif
5527       enddo
5528       if (lprn) then
5529         write (iout,*) "stochforcvec"
5530         do i=1,3*dimen
5531           write(iout,'(i5,e15.5)') i,stochforcvec(i)
5532         enddo
5533         write(iout,*) "Stochastic forces backbone"
5534         do i=0,nct-1
5535           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5536         enddo
5537         write(iout,*) "Stochastic forces side chain"
5538         do i=nnt,nct
5539           write(iout,'(i5,3e15.5)') &
5540             i,(stochforc(j,i+nres),j=1,3)
5541         enddo
5542       endif
5543
5544       if (lprn) then
5545
5546       ind=0
5547       do i=nnt,nct-1
5548         write (iout,*) i,ind
5549         do j=1,3
5550           forcvec(ind+j)=force(j,i)
5551         enddo
5552         ind=ind+3
5553       enddo
5554       do i=nnt,nct
5555         write (iout,*) i,ind
5556         do j=1,3
5557           forcvec(j+ind)=force(j,i+nres)
5558         enddo
5559         ind=ind+3
5560       enddo
5561
5562       write (iout,*) "forcvec"
5563       ind=0
5564       do i=nnt,nct-1
5565         do j=1,3
5566           write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5567             forcvec(ind+j)
5568         enddo
5569         ind=ind+3
5570       enddo
5571       do i=nnt,nct
5572         do j=1,3
5573           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5574            forcvec(ind+j)
5575         enddo
5576         ind=ind+3
5577       enddo
5578
5579       endif
5580
5581       return
5582       end subroutine stochastic_force
5583 !-----------------------------------------------------------------------------
5584       subroutine sdarea(gamvec)
5585 !
5586 ! Scale the friction coefficients according to solvent accessible surface areas
5587 ! Code adapted from TINKER
5588 ! AL 9/3/04
5589 !
5590       use energy_data
5591 !      implicit real*8 (a-h,o-z)
5592 !      include 'DIMENSIONS'
5593 !      include 'COMMON.CONTROL'
5594 !      include 'COMMON.VAR'
5595 !      include 'COMMON.MD'
5596 !#ifndef LANG0
5597 !      include 'COMMON.LANGEVIN'
5598 !#else
5599 !      include 'COMMON.LANGEVIN.lang0'
5600 !#endif
5601 !      include 'COMMON.CHAIN'
5602 !      include 'COMMON.DERIV'
5603 !      include 'COMMON.GEO'
5604 !      include 'COMMON.LOCAL'
5605 !      include 'COMMON.INTERACT'
5606 !      include 'COMMON.IOUNITS'
5607 !      include 'COMMON.NAMES'
5608       real(kind=8),dimension(2*nres) :: radius,gamvec   !(maxres2)
5609       real(kind=8),parameter :: twosix = 1.122462048309372981d0
5610       logical :: lprn = .false.
5611       real(kind=8) :: probe,area,ratio
5612       integer :: i,j,ind,iti,mnum
5613 !
5614 !     determine new friction coefficients every few SD steps
5615 !
5616 !     set the atomic radii to estimates of sigma values
5617 !
5618 !      print *,"Entered sdarea"
5619       probe = 0.0d0
5620
5621       do i=1,2*nres
5622         radius(i)=0.0d0
5623       enddo
5624 !  Load peptide group radii
5625       do i=nnt,nct-1
5626         mnum=molnum(i)
5627         radius(i)=pstok(mnum)
5628       enddo
5629 !  Load side chain radii
5630       do i=nnt,nct
5631         mnum=molnum(i)
5632         iti=itype(i,mnum)
5633         radius(i+nres)=restok(iti,mnum)
5634       enddo
5635 !      do i=1,2*nres
5636 !        write (iout,*) "i",i," radius",radius(i)
5637 !      enddo
5638       do i = 1, 2*nres
5639          radius(i) = radius(i) / twosix
5640          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
5641       end do
5642 !
5643 !     scale atomic friction coefficients by accessible area
5644 !
5645       if (lprn) write (iout,*) &
5646         "Original gammas, surface areas, scaling factors, new gammas, ",&
5647         "std's of stochastic forces"
5648       ind=0
5649       do i=nnt,nct-1
5650         if (radius(i).gt.0.0d0) then
5651           call surfatom (i,area,radius)
5652           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5653           if (lprn) write (iout,'(i5,3f10.5,$)') &
5654             i,gamvec(ind+1),area,ratio
5655           do j=1,3
5656             ind=ind+1
5657             gamvec(ind) = ratio * gamvec(ind)
5658           enddo
5659           mnum=molnum(i)
5660           stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5661           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5662         endif
5663       enddo
5664       do i=nnt,nct
5665         if (radius(i+nres).gt.0.0d0) then
5666           call surfatom (i+nres,area,radius)
5667           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5668           if (lprn) write (iout,'(i5,3f10.5,$)') &
5669             i,gamvec(ind+1),area,ratio
5670           do j=1,3
5671             ind=ind+1
5672             gamvec(ind) = ratio * gamvec(ind)
5673           enddo
5674           mnum=molnum(i)
5675           stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5676           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5677         endif
5678       enddo
5679
5680       return
5681       end subroutine sdarea
5682 !-----------------------------------------------------------------------------
5683 ! surfatom.f
5684 !-----------------------------------------------------------------------------
5685 !
5686 !
5687 !     ###################################################
5688 !     ##  COPYRIGHT (C)  1996  by  Jay William Ponder  ##
5689 !     ##              All Rights Reserved              ##
5690 !     ###################################################
5691 !
5692 !     ################################################################
5693 !     ##                                                            ##
5694 !     ##  subroutine surfatom  --  exposed surface area of an atom  ##
5695 !     ##                                                            ##
5696 !     ################################################################
5697 !
5698 !
5699 !     "surfatom" performs an analytical computation of the surface
5700 !     area of a specified atom; a simplified version of "surface"
5701 !
5702 !     literature references:
5703 !
5704 !     T. J. Richmond, "Solvent Accessible Surface Area and
5705 !     Excluded Volume in Proteins", Journal of Molecular Biology,
5706 !     178, 63-89 (1984)
5707 !
5708 !     L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5709 !     Applied to Molecular Dynamics of Proteins in Solution",
5710 !     Protein Science, 1, 227-235 (1992)
5711 !
5712 !     variables and parameters:
5713 !
5714 !     ir       number of atom for which area is desired
5715 !     area     accessible surface area of the atom
5716 !     radius   radii of each of the individual atoms
5717 !
5718 !
5719       subroutine surfatom(ir,area,radius)
5720
5721 !      implicit real*8 (a-h,o-z)
5722 !      include 'DIMENSIONS'
5723 !      include 'sizes.i'
5724 !      include 'COMMON.GEO'
5725 !      include 'COMMON.IOUNITS'
5726 !      integer :: nres,
5727       integer :: nsup,nstart_sup
5728 !      double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5729 !      common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5730 !     & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5731 !     & dc_work(MAXRES6),nres,nres0
5732       integer,parameter :: maxarc=300
5733       integer :: i,j,k,m
5734       integer :: ii,ib,jb
5735       integer :: io,ir
5736       integer :: mi,ni,narc
5737       integer :: key(maxarc)
5738       integer :: intag(maxarc)
5739       integer :: intag1(maxarc)
5740       real(kind=8) :: area,arcsum
5741       real(kind=8) :: arclen,exang
5742       real(kind=8) :: delta,delta2
5743       real(kind=8) :: eps,rmove
5744       real(kind=8) :: xr,yr,zr
5745       real(kind=8) :: rr,rrsq
5746       real(kind=8) :: rplus,rminus
5747       real(kind=8) :: axx,axy,axz
5748       real(kind=8) :: ayx,ayy
5749       real(kind=8) :: azx,azy,azz
5750       real(kind=8) :: uxj,uyj,uzj
5751       real(kind=8) :: tx,ty,tz
5752       real(kind=8) :: txb,tyb,td
5753       real(kind=8) :: tr2,tr,txr,tyr
5754       real(kind=8) :: tk1,tk2
5755       real(kind=8) :: thec,the,t,tb
5756       real(kind=8) :: txk,tyk,tzk
5757       real(kind=8) :: t1,ti,tf,tt
5758       real(kind=8) :: txj,tyj,tzj
5759       real(kind=8) :: ccsq,cc,xysq
5760       real(kind=8) :: bsqk,bk,cosine
5761       real(kind=8) :: dsqj,gi,pix2
5762       real(kind=8) :: therk,dk,gk
5763       real(kind=8) :: risqk,rik
5764       real(kind=8) :: radius(2*nres)    !(maxatm) (maxatm=maxres2)
5765       real(kind=8) :: ri(maxarc),risq(maxarc)
5766       real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5767       real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5768       real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5769       real(kind=8) :: dsq(maxarc),bsq(maxarc)
5770       real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5771       real(kind=8) :: arci(maxarc),arcf(maxarc)
5772       real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5773       real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5774       real(kind=8) :: kent(maxarc),kout(maxarc)
5775       real(kind=8) :: ther(maxarc)
5776       logical :: moved,top
5777       logical :: omit(maxarc)
5778 !
5779 !      include 'sizes.i'
5780 !      maxatm = 2*nres  !maxres2 maxres2=2*maxres
5781 !      maxlight = 8*maxatm
5782 !      maxbnd = 2*maxatm
5783 !      maxang = 3*maxatm
5784 !      maxtors = 4*maxatm
5785 !
5786 !     zero out the surface area for the sphere of interest
5787 !
5788       area = 0.0d0
5789 !      write (2,*) "ir",ir," radius",radius(ir)
5790       if (radius(ir) .eq. 0.0d0)  return
5791 !
5792 !     set the overlap significance and connectivity shift
5793 !
5794       pix2 = 2.0d0 * pi
5795       delta = 1.0d-8
5796       delta2 = delta * delta
5797       eps = 1.0d-8
5798       moved = .false.
5799       rmove = 1.0d-8
5800 !
5801 !     store coordinates and radius of the sphere of interest
5802 !
5803       xr = c(1,ir)
5804       yr = c(2,ir)
5805       zr = c(3,ir)
5806       rr = radius(ir)
5807       rrsq = rr * rr
5808 !
5809 !     initialize values of some counters and summations
5810 !
5811    10 continue
5812       io = 0
5813       jb = 0
5814       ib = 0
5815       arclen = 0.0d0
5816       exang = 0.0d0
5817 !
5818 !     test each sphere to see if it overlaps the sphere of interest
5819 !
5820       do i = 1, 2*nres
5821          if (i.eq.ir .or. radius(i).eq.0.0d0)  goto 30
5822          rplus = rr + radius(i)
5823          tx = c(1,i) - xr
5824          if (abs(tx) .ge. rplus)  goto 30
5825          ty = c(2,i) - yr
5826          if (abs(ty) .ge. rplus)  goto 30
5827          tz = c(3,i) - zr
5828          if (abs(tz) .ge. rplus)  goto 30
5829 !
5830 !     check for sphere overlap by testing distance against radii
5831 !
5832          xysq = tx*tx + ty*ty
5833          if (xysq .lt. delta2) then
5834             tx = delta
5835             ty = 0.0d0
5836             xysq = delta2
5837          end if
5838          ccsq = xysq + tz*tz
5839          cc = sqrt(ccsq)
5840          if (rplus-cc .le. delta)  goto 30
5841          rminus = rr - radius(i)
5842 !
5843 !     check to see if sphere of interest is completely buried
5844 !
5845          if (cc-abs(rminus) .le. delta) then
5846             if (rminus .le. 0.0d0)  goto 170
5847             goto 30
5848          end if
5849 !
5850 !     check for too many overlaps with sphere of interest
5851 !
5852          if (io .ge. maxarc) then
5853             write (iout,20)
5854    20       format (/,' SURFATOM  --  Increase the Value of MAXARC')
5855             stop
5856          end if
5857 !
5858 !     get overlap between current sphere and sphere of interest
5859 !
5860          io = io + 1
5861          xc1(io) = tx
5862          yc1(io) = ty
5863          zc1(io) = tz
5864          dsq1(io) = xysq
5865          bsq1(io) = ccsq
5866          b1(io) = cc
5867          gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5868          intag1(io) = i
5869          omit(io) = .false.
5870    30    continue
5871       end do
5872 !
5873 !     case where no other spheres overlap the sphere of interest
5874 !
5875       if (io .eq. 0) then
5876          area = 4.0d0 * pi * rrsq
5877          return
5878       end if
5879 !
5880 !     case where only one sphere overlaps the sphere of interest
5881 !
5882       if (io .eq. 1) then
5883          area = pix2 * (1.0d0 + gr(1))
5884          area = mod(area,4.0d0*pi) * rrsq
5885          return
5886       end if
5887 !
5888 !     case where many spheres intersect the sphere of interest;
5889 !     sort the intersecting spheres by their degree of overlap
5890 !
5891       call sort2 (io,gr,key)
5892       do i = 1, io
5893          k = key(i)
5894          intag(i) = intag1(k)
5895          xc(i) = xc1(k)
5896          yc(i) = yc1(k)
5897          zc(i) = zc1(k)
5898          dsq(i) = dsq1(k)
5899          b(i) = b1(k)
5900          bsq(i) = bsq1(k)
5901       end do
5902 !
5903 !     get radius of each overlap circle on surface of the sphere
5904 !
5905       do i = 1, io
5906          gi = gr(i) * rr
5907          bg(i) = b(i) * gi
5908          risq(i) = rrsq - gi*gi
5909          ri(i) = sqrt(risq(i))
5910          ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5911       end do
5912 !
5913 !     find boundary of inaccessible area on sphere of interest
5914 !
5915       do k = 1, io-1
5916          if (.not. omit(k)) then
5917             txk = xc(k)
5918             tyk = yc(k)
5919             tzk = zc(k)
5920             bk = b(k)
5921             therk = ther(k)
5922 !
5923 !     check to see if J circle is intersecting K circle;
5924 !     get distance between circle centers and sum of radii
5925 !
5926             do j = k+1, io
5927                if (omit(j))  goto 60
5928                cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5929                cc = acos(min(1.0d0,max(-1.0d0,cc)))
5930                td = therk + ther(j)
5931 !
5932 !     check to see if circles enclose separate regions
5933 !
5934                if (cc .ge. td)  goto 60
5935 !
5936 !     check for circle J completely inside circle K
5937 !
5938                if (cc+ther(j) .lt. therk)  goto 40
5939 !
5940 !     check for circles that are essentially parallel
5941 !
5942                if (cc .gt. delta)  goto 50
5943    40          continue
5944                omit(j) = .true.
5945                goto 60
5946 !
5947 !     check to see if sphere of interest is completely buried
5948 !
5949    50          continue
5950                if (pix2-cc .le. td)  goto 170
5951    60          continue
5952             end do
5953          end if
5954       end do
5955 !
5956 !     find T value of circle intersections
5957 !
5958       do k = 1, io
5959          if (omit(k))  goto 110
5960          omit(k) = .true.
5961          narc = 0
5962          top = .false.
5963          txk = xc(k)
5964          tyk = yc(k)
5965          tzk = zc(k)
5966          dk = sqrt(dsq(k))
5967          bsqk = bsq(k)
5968          bk = b(k)
5969          gk = gr(k) * rr
5970          risqk = risq(k)
5971          rik = ri(k)
5972          therk = ther(k)
5973 !
5974 !     rotation matrix elements
5975 !
5976          t1 = tzk / (bk*dk)
5977          axx = txk * t1
5978          axy = tyk * t1
5979          axz = dk / bk
5980          ayx = tyk / dk
5981          ayy = txk / dk
5982          azx = txk / bk
5983          azy = tyk / bk
5984          azz = tzk / bk
5985          do j = 1, io
5986             if (.not. omit(j)) then
5987                txj = xc(j)
5988                tyj = yc(j)
5989                tzj = zc(j)
5990 !
5991 !     rotate spheres so K vector colinear with z-axis
5992 !
5993                uxj = txj*axx + tyj*axy - tzj*axz
5994                uyj = tyj*ayy - txj*ayx
5995                uzj = txj*azx + tyj*azy + tzj*azz
5996                cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5997                if (acos(cosine) .lt. therk+ther(j)) then
5998                   dsqj = uxj*uxj + uyj*uyj
5999                   tb = uzj*gk - bg(j)
6000                   txb = uxj * tb
6001                   tyb = uyj * tb
6002                   td = rik * dsqj
6003                   tr2 = risqk*dsqj - tb*tb
6004                   tr2 = max(eps,tr2)
6005                   tr = sqrt(tr2)
6006                   txr = uxj * tr
6007                   tyr = uyj * tr
6008 !
6009 !     get T values of intersection for K circle
6010 !
6011                   tb = (txb+tyr) / td
6012                   tb = min(1.0d0,max(-1.0d0,tb))
6013                   tk1 = acos(tb)
6014                   if (tyb-txr .lt. 0.0d0)  tk1 = pix2 - tk1
6015                   tb = (txb-tyr) / td
6016                   tb = min(1.0d0,max(-1.0d0,tb))
6017                   tk2 = acos(tb)
6018                   if (tyb+txr .lt. 0.0d0)  tk2 = pix2 - tk2
6019                   thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6020                   if (abs(thec) .lt. 1.0d0) then
6021                      the = -acos(thec)
6022                   else if (thec .ge. 1.0d0) then
6023                      the = 0.0d0
6024                   else if (thec .le. -1.0d0) then
6025                      the = -pi
6026                   end if
6027 !
6028 !     see if "tk1" is entry or exit point; check t=0 point;
6029 !     "ti" is exit point, "tf" is entry point
6030 !
6031                   cosine = min(1.0d0,max(-1.0d0, &
6032                                   (uzj*gk-uxj*rik)/(b(j)*rr)))
6033                   if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6034                      ti = tk2
6035                      tf = tk1
6036                   else
6037                      ti = tk2
6038                      tf = tk1
6039                   end if
6040                   narc = narc + 1
6041                   if (narc .ge. maxarc) then
6042                      write (iout,70)
6043    70                format (/,' SURFATOM  --  Increase the Value',&
6044                                 ' of MAXARC')
6045                      stop
6046                   end if
6047                   if (tf .le. ti) then
6048                      arcf(narc) = tf
6049                      arci(narc) = 0.0d0
6050                      tf = pix2
6051                      lt(narc) = j
6052                      ex(narc) = the
6053                      top = .true.
6054                      narc = narc + 1
6055                   end if
6056                   arcf(narc) = tf
6057                   arci(narc) = ti
6058                   lt(narc) = j
6059                   ex(narc) = the
6060                   ux(j) = uxj
6061                   uy(j) = uyj
6062                   uz(j) = uzj
6063                end if
6064             end if
6065          end do
6066          omit(k) = .false.
6067 !
6068 !     special case; K circle without intersections
6069 !
6070          if (narc .le. 0)  goto 90
6071 !
6072 !     general case; sum up arclength and set connectivity code
6073 !
6074          call sort2 (narc,arci,key)
6075          arcsum = arci(1)
6076          mi = key(1)
6077          t = arcf(mi)
6078          ni = mi
6079          if (narc .gt. 1) then
6080             do j = 2, narc
6081                m = key(j)
6082                if (t .lt. arci(j)) then
6083                   arcsum = arcsum + arci(j) - t
6084                   exang = exang + ex(ni)
6085                   jb = jb + 1
6086                   if (jb .ge. maxarc) then
6087                      write (iout,80)
6088    80                format (/,' SURFATOM  --  Increase the Value',&
6089                                 ' of MAXARC')
6090                      stop
6091                   end if
6092                   i = lt(ni)
6093                   kent(jb) = maxarc*i + k
6094                   i = lt(m)
6095                   kout(jb) = maxarc*k + i
6096                end if
6097                tt = arcf(m)
6098                if (tt .ge. t) then
6099                   t = tt
6100                   ni = m
6101                end if
6102             end do
6103          end if
6104          arcsum = arcsum + pix2 - t
6105          if (.not. top) then
6106             exang = exang + ex(ni)
6107             jb = jb + 1
6108             i = lt(ni)
6109             kent(jb) = maxarc*i + k
6110             i = lt(mi)
6111             kout(jb) = maxarc*k + i
6112          end if
6113          goto 100
6114    90    continue
6115          arcsum = pix2
6116          ib = ib + 1
6117   100    continue
6118          arclen = arclen + gr(k)*arcsum
6119   110    continue
6120       end do
6121       if (arclen .eq. 0.0d0)  goto 170
6122       if (jb .eq. 0)  goto 150
6123 !
6124 !     find number of independent boundaries and check connectivity
6125 !
6126       j = 0
6127       do k = 1, jb
6128          if (kout(k) .ne. 0) then
6129             i = k
6130   120       continue
6131             m = kout(i)
6132             kout(i) = 0
6133             j = j + 1
6134             do ii = 1, jb
6135                if (m .eq. kent(ii)) then
6136                   if (ii .eq. k) then
6137                      ib = ib + 1
6138                      if (j .eq. jb)  goto 150
6139                      goto 130
6140                   end if
6141                   i = ii
6142                   goto 120
6143                end if
6144             end do
6145   130       continue
6146          end if
6147       end do
6148       ib = ib + 1
6149 !
6150 !     attempt to fix connectivity error by moving atom slightly
6151 !
6152       if (moved) then
6153          write (iout,140)  ir
6154   140    format (/,' SURFATOM  --  Connectivity Error at Atom',i6)
6155       else
6156          moved = .true.
6157          xr = xr + rmove
6158          yr = yr + rmove
6159          zr = zr + rmove
6160          goto 10
6161       end if
6162 !
6163 !     compute the exposed surface area for the sphere of interest
6164 !
6165   150 continue
6166       area = ib*pix2 + exang + arclen
6167       area = mod(area,4.0d0*pi) * rrsq
6168 !
6169 !     attempt to fix negative area by moving atom slightly
6170 !
6171       if (area .lt. 0.0d0) then
6172          if (moved) then
6173             write (iout,160)  ir
6174   160       format (/,' SURFATOM  --  Negative Area at Atom',i6)
6175          else
6176             moved = .true.
6177             xr = xr + rmove
6178             yr = yr + rmove
6179             zr = zr + rmove
6180             goto 10
6181          end if
6182       end if
6183   170 continue
6184       return
6185       end subroutine surfatom
6186 !----------------------------------------------------------------
6187 !----------------------------------------------------------------
6188       subroutine alloc_MD_arrays
6189 !EL Allocation of arrays used by MD module
6190
6191       integer :: nres2,nres6
6192       nres2=nres*2
6193       nres6=nres*6
6194 !----------------------
6195 #ifndef LANG0
6196 ! commom.langevin
6197 !      common /langforc/
6198       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6199       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6200       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6201       allocate(fricvec(nres2,nres2))
6202       allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6203       allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6204       allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6205       allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6206       allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6207       allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6208       allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6209       allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6210       allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6211       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6212 !      common /langmat/
6213       allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6214 !----------------------
6215 #else
6216 ! commom.langevin.lang0
6217 !      common /langforc/
6218       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6219       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6220       allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6221       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6222       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6223 #endif
6224
6225       if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6226 !----------------------
6227 ! commom.hairpin in CSA module
6228 !----------------------
6229 ! common.mce in MCM_MD module
6230 !----------------------
6231 ! common.MD
6232 !      common /mdgrad/ in module.energy
6233 !      common /back_constr/ in module.energy
6234 !      common /qmeas/ in module.energy
6235 !      common /mdpar/
6236 !      common /MDcalc/
6237       allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6238 !      common /lagrange/
6239       allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6240       allocate(d_a_work(nres6)) !(6*MAXRES)
6241 #ifdef FIVEDIAG
6242       allocate(DM(nres2),DU1(nres2),DU2(nres2))
6243       allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6244 #else
6245       write (iout,*) "Before A Allocation",nfgtasks-1
6246       call flush(iout)
6247       allocate(Gmat(nres2,nres2),A(nres2,nres2))
6248       if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6249       allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6250 #endif
6251       allocate(Geigen(nres2)) !(maxres2)
6252       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6253 !      common /inertia/ in io_conf: parmread
6254 !      real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6255 !      common /langevin/in io read_MDpar
6256 !      real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6257 !      real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6258 ! in io_conf: parmread
6259 !      real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6260 !      common /mdpmpi/ in control: ergastulum
6261       if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6262       if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6263       if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6264       if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6265 !----------------------
6266 ! common.muca in read_muca
6267 !      common /double_muca/
6268 !      real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6269 !      real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6270 !       nemuca2,hist !(4*maxres)
6271 !      common /integer_muca/
6272 !      integer :: nmuca,imtime,muca_smooth
6273 !      common /mucarem/
6274 !      real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6275 !----------------------
6276 ! common.MD
6277 !      common /mdgrad/ in module.energy
6278 !      common /back_constr/ in module.energy
6279 !      common /qmeas/ in module.energy
6280 !      common /mdpar/
6281 !      common /MDcalc/
6282 !      common /lagrange/
6283       allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6284       allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6285       allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6286       allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6287 !----------------------
6288 !      COMMON /BANII/ D
6289       allocate(D_ban(nres6))    !(MAXRES6) maxres6=6*maxres
6290 !      common /stochcalc/ stochforcvec
6291       allocate(stochforcvec(nres6))     !(MAXRES6) maxres6=6*maxres
6292 !----------------------
6293       return
6294       end subroutine alloc_MD_arrays
6295 !-----------------------------------------------------------------------------
6296 !-----------------------------------------------------------------------------
6297       end module MDyn