source/unres/MD.F90

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