source/unres/MD.f90

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