source/unres/MD.F90

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