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