source/unres/MD.f90

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