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