source/unres/MD.F90

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