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