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.99e18 .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) then
2511 !         call chainbuild
2512          if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2513           if (overlapsc) then
2514            print *, 'Calling OVERLAP_SC'
2515            call overlap_sc(fail)
2516            print *,'after OVERLAP'
2517           endif
2518           if (searchsc) then
2519            print *,'call SC_MOVE'
2520            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2521            print *,'SC_move',nft_sc,etot
2522            if(me.eq.king.or..not.out1file) &
2523             write(iout,*) 'SC_move',nft_sc,etot
2524           endif
2525
2526           if(dccart)then
2527            print *, 'Calling MINIM_DC'
2528            call minim_dc(etot,iretcode,nfun)
2529           else
2530            call geom_to_var(nvar,varia)
2531            print *,'Calling MINIMIZE.'
2532            call minimize(etot,varia,iretcode,nfun)
2533            call var_to_geom(nvar,varia)
2534           endif
2535           if(me.eq.king.or..not.out1file) &
2536              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2537          endif
2538       endif
2539       call chainbuild_cart
2540       call kinetic(EK)
2541       if (tbf) then
2542         call verlet_bath
2543       endif
2544       kinetic_T=2.0d0/(dimen3*Rb)*EK
2545       if(me.eq.king.or..not.out1file)then
2546        call cartprint
2547        call intout
2548       endif
2549 #ifdef MPI
2550       tt0=MPI_Wtime()
2551 #else
2552       tt0=tcpu()
2553 #endif
2554       call zerograd
2555       write(iout,*) "before ETOTAL"
2556       call etotal(potEcomp)
2557       if (large) call enerprint(potEcomp)
2558 #ifdef TIMING_ENE
2559 #ifdef MPI
2560       t_etotal=t_etotal+MPI_Wtime()-tt0
2561 #else
2562       t_etotal=t_etotal+tcpu()-tt0
2563 #endif
2564 #endif
2565       potE=potEcomp(0)
2566       call cartgrad
2567       call lagrangian
2568       call max_accel
2569       if (amax*d_time .gt. dvmax) then
2570         d_time=d_time*dvmax/amax
2571         if(me.eq.king.or..not.out1file) write (iout,*) &
2572          "Time step reduced to",d_time,&
2573          " because of too large initial acceleration."
2574       endif
2575       if(me.eq.king.or..not.out1file)then
2576        write(iout,*) "Potential energy and its components"
2577        call enerprint(potEcomp)
2578 !       write(iout,*) (potEcomp(i),i=0,n_ene)
2579       endif
2580       potE=potEcomp(0)-potEcomp(20)
2581       totE=EK+potE
2582       itime=0
2583       if (ntwe.ne.0) call statout(itime)
2584       if(me.eq.king.or..not.out1file) &
2585         write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2586          " Kinetic energy",EK," Potential energy",potE, &
2587          " Total energy",totE," Maximum acceleration ", &
2588          amax
2589       if (large) then
2590         write (iout,*) "Initial coordinates"
2591         do i=1,nres
2592           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2593           (c(j,i+nres),j=1,3)
2594         enddo
2595         write (iout,*) "Initial dC"
2596         do i=0,nres
2597           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2598           (dc(j,i+nres),j=1,3)
2599         enddo
2600         write (iout,*) "Initial velocities"
2601         write (iout,"(13x,' backbone ',23x,' side chain')")
2602         do i=0,nres
2603           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2604           (d_t(j,i+nres),j=1,3)
2605         enddo
2606         write (iout,*) "Initial accelerations"
2607         do i=0,nres
2608 !          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2609           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2610           (d_a(j,i+nres),j=1,3)
2611         enddo
2612       endif
2613       do i=0,2*nres
2614         do j=1,3
2615           dc_old(j,i)=dc(j,i)
2616           d_t_old(j,i)=d_t(j,i)
2617           d_a_old(j,i)=d_a(j,i)
2618         enddo
2619 !        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2620       enddo
2621       if (RESPA) then
2622 #ifdef MPI
2623         tt0 =MPI_Wtime()
2624 #else
2625         tt0 = tcpu()
2626 #endif
2627         call zerograd
2628         call etotal_short(energia_short)
2629         if (large) call enerprint(potEcomp)
2630 #ifdef TIMING_ENE
2631 #ifdef MPI
2632         t_eshort=t_eshort+MPI_Wtime()-tt0
2633 #else
2634         t_eshort=t_eshort+tcpu()-tt0
2635 #endif
2636 #endif
2637         call cartgrad
2638         call lagrangian
2639         if(.not.out1file .and. large) then
2640           write (iout,*) "energia_long",energia_long(0),&
2641            " energia_short",energia_short(0),&
2642            " total",energia_long(0)+energia_short(0)
2643           write (iout,*) "Initial fast-force accelerations"
2644           do i=0,nres
2645             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2646             (d_a(j,i+nres),j=1,3)
2647           enddo
2648         endif
2649 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2650         do i=0,2*nres
2651           do j=1,3
2652             d_a_short(j,i)=d_a(j,i)
2653           enddo
2654         enddo
2655 #ifdef MPI
2656         tt0=MPI_Wtime()
2657 #else
2658         tt0=tcpu()
2659 #endif
2660         call zerograd
2661         call etotal_long(energia_long)
2662         if (large) call enerprint(potEcomp)
2663 #ifdef TIMING_ENE
2664 #ifdef MPI
2665         t_elong=t_elong+MPI_Wtime()-tt0
2666 #else
2667         t_elong=t_elong+tcpu()-tt0
2668 #endif
2669 #endif
2670         call cartgrad
2671         call lagrangian
2672         if(.not.out1file .and. large) then
2673           write (iout,*) "energia_long",energia_long(0)
2674           write (iout,*) "Initial slow-force accelerations"
2675           do i=0,nres
2676             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2677             (d_a(j,i+nres),j=1,3)
2678           enddo
2679         endif
2680 #ifdef MPI
2681         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2682 #else
2683         t_enegrad=t_enegrad+tcpu()-tt0
2684 #endif
2685       endif
2686       return
2687       end subroutine init_MD
2688 !-----------------------------------------------------------------------------
2689       subroutine random_vel
2690
2691 !      implicit real*8 (a-h,o-z)
2692       use energy_data
2693       use random, only:anorm_distr
2694       use MD_data
2695 !      include 'DIMENSIONS'
2696 !      include 'COMMON.CONTROL'
2697 !      include 'COMMON.VAR'
2698 !      include 'COMMON.MD'
2699 !#ifndef LANG0
2700 !      include 'COMMON.LANGEVIN'
2701 !#else
2702 !      include 'COMMON.LANGEVIN.lang0'
2703 !#endif
2704 !      include 'COMMON.CHAIN'
2705 !      include 'COMMON.DERIV'
2706 !      include 'COMMON.GEO'
2707 !      include 'COMMON.LOCAL'
2708 !      include 'COMMON.INTERACT'
2709 !      include 'COMMON.IOUNITS'
2710 !      include 'COMMON.NAMES'
2711 !      include 'COMMON.TIME1'
2712       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
2713 !#define DEBUG
2714 #ifdef FIVEDIAG
2715        real(kind=8) ,allocatable, dimension(:)  :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2716        real(kind=8) :: sumx
2717 #ifdef DEBUG
2718        real(kind=8) ,allocatable, dimension(:)  :: rsold
2719        real (kind=8),allocatable,dimension(:,:) :: matold
2720        integer :: iti
2721 #endif
2722 #endif
2723       integer :: i,j,ii,k,ind,mark,imark,mnum
2724 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2725 ! First generate velocities in the eigenspace of the G matrix
2726 !      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2727 !      call flush(iout)
2728       xv=0.0d0
2729       ii=0
2730       do i=1,dimen
2731         do k=1,3
2732           ii=ii+1
2733           sigv=dsqrt((Rb*t_bath)/geigen(i))
2734           lowb=-5*sigv
2735           highb=5*sigv
2736           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2737 #ifdef DEBUG
2738           write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2739             " d_t_work_new",d_t_work_new(ii)
2740 #endif
2741         enddo
2742       enddo
2743 #ifdef DEBUG
2744 ! diagnostics
2745       Ek1=0.0d0
2746       ii=0
2747       do i=1,dimen
2748         do k=1,3
2749           ii=ii+1
2750           Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2751         enddo
2752       enddo
2753       write (iout,*) "Ek from eigenvectors",Ek1
2754       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2755 ! end diagnostics
2756 #endif
2757 #ifdef FIVEDIAG
2758 ! Transform velocities to UNRES coordinate space
2759        allocate (DL1(2*nres))
2760        allocate (DDU1(2*nres))
2761        allocate (DL2(2*nres))
2762        allocate (DDU2(2*nres))
2763        allocate (xsolv(2*nres))
2764        allocate (DML(2*nres))
2765        allocate (rs(2*nres))
2766 #ifdef DEBUG
2767        allocate (rsold(2*nres))
2768        allocate (matold(2*nres,2*nres))
2769        matold=0
2770        matold(1,1)=DMorig(1)
2771        matold(1,2)=DU1orig(1)
2772        matold(1,3)=DU2orig(1)
2773        write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2774
2775       do i=2,dimen-1
2776         do j=1,dimen
2777           if (i.eq.j) then
2778              matold(i,j)=DMorig(i)
2779              matold(i,j-1)=DU1orig(i-1)
2780            if (j.ge.3) then
2781                  matold(i,j-2)=DU2orig(i-2)
2782                endif
2783
2784             endif
2785
2786           enddo
2787           do j=1,dimen-1
2788             if (i.eq.j) then
2789               matold(i,j+1)=DU1orig(i)
2790
2791              end if
2792           enddo
2793           do j=1,dimen-2
2794             if(i.eq.j) then
2795                matold(i,j+2)=DU2orig(i)
2796             endif
2797           enddo
2798        enddo
2799        matold(dimen,dimen)=DMorig(dimen)
2800        matold(dimen,dimen-1)=DU1orig(dimen-1)
2801        matold(dimen,dimen-2)=DU2orig(dimen-2)
2802        write(iout,*) "old gmatrix"
2803        call matout(dimen,dimen,2*nres,2*nres,matold)
2804 #endif
2805       d_t_work=0.0d0
2806       do i=1,dimen
2807 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2808        do imark=dimen,1,-1
2809
2810          do j=1,imark-1
2811            DML(j)=DMorig(j)-geigen(i)
2812          enddo
2813          do j=imark+1,dimen
2814            DML(j-1)=DMorig(j)-geigen(i)
2815          enddo
2816          do j=1,imark-2
2817            DDU1(j)=DU1orig(j)
2818          enddo
2819          DDU1(imark-1)=DU2orig(imark-1)
2820          do j=imark+1,dimen-1
2821            DDU1(j-1)=DU1orig(j)
2822          enddo
2823          do j=1,imark-3
2824            DDU2(j)=DU2orig(j)
2825          enddo
2826          DDU2(imark-2)=0.0d0
2827          DDU2(imark-1)=0.0d0
2828          do j=imark,dimen-3
2829            DDU2(j)=DU2orig(j+1)
2830          enddo
2831          do j=1,dimen-3
2832            DL2(j+2)=DDU2(j)
2833          enddo
2834          do j=1,dimen-2
2835            DL1(j+1)=DDU1(j)
2836          enddo
2837 #ifdef DEBUG
2838          write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2839          write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2840          write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2841          write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2842          write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2843          write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2844 #endif
2845          rs=0.0d0
2846          if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2847          if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2848          if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2849          if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2850 #ifdef DEBUG
2851          rsold=rs
2852 #endif
2853 !         write (iout,*) "Vector rs"
2854 !         do j=1,dimen-1
2855 !           write (iout,*) j,rs(j)
2856 !         enddo
2857          xsolv=0.0d0
2858          call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2859
2860          if (mark.eq.1) then
2861
2862 #ifdef DEBUG
2863 ! Check solution
2864          do j=1,imark-1
2865            sumx=-geigen(i)*xsolv(j)
2866            do k=1,imark-1
2867              sumx=sumx+matold(j,k)*xsolv(k)
2868            enddo
2869            do k=imark+1,dimen
2870              sumx=sumx+matold(j,k)*xsolv(k-1)
2871            enddo
2872            write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2873          enddo
2874          do j=imark+1,dimen
2875            sumx=-geigen(i)*xsolv(j-1)
2876            do k=1,imark-1
2877              sumx=sumx+matold(j,k)*xsolv(k)
2878            enddo
2879            do k=imark+1,dimen
2880              sumx=sumx+matold(j,k)*xsolv(k-1)
2881            enddo
2882            write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2883          enddo
2884 ! end check
2885          write (iout,*)&
2886           "Solution of equations system",i," for eigenvalue",geigen(i)
2887          do j=1,dimen-1
2888            write(iout,'(i5,f10.5)') j,xsolv(j)
2889          enddo
2890 #endif
2891          do j=dimen-1,imark,-1
2892            xsolv(j+1)=xsolv(j)
2893          enddo
2894          xsolv(imark)=1.0d0
2895 #ifdef DEBUG
2896          write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2897          do j=1,dimen
2898            write(iout,'(i5,f10.5)') j,xsolv(j)
2899          enddo
2900 #endif
2901 ! Normalize ith eigenvector
2902          sumx=0.0d0
2903          do j=1,dimen
2904            sumx=sumx+xsolv(j)**2
2905          enddo
2906          sumx=dsqrt(sumx)
2907          do j=1,dimen
2908            xsolv(j)=xsolv(j)/sumx
2909          enddo
2910 #ifdef DEBUG
2911          write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2912          do j=1,dimen
2913            write(iout,'(i5,3f10.5)') j,xsolv(j)
2914          enddo
2915 #endif
2916 ! All done at this point for eigenvector i, exit loop
2917          exit
2918
2919          endif ! mark.eq.1
2920
2921        enddo ! imark
2922
2923        if (mark.ne.1) then
2924          write (iout,*) "Unable to find eigenvector",i
2925        endif
2926
2927 !       write (iout,*) "i=",i
2928        do k=1,3
2929 !         write (iout,*) "k=",k
2930          do j=1,dimen
2931            ind=(j-1)*3+k
2932 !           write(iout,*) "ind",ind," ind1",3*(i-1)+k
2933            d_t_work(ind)=d_t_work(ind) &
2934                             +xsolv(j)*d_t_work_new(3*(i-1)+k)
2935          enddo
2936        enddo
2937       enddo ! i (loop over the eigenvectors)
2938
2939 #ifdef DEBUG
2940       write (iout,*) "d_t_work"
2941       do i=1,3*dimen
2942         write (iout,"(i5,f10.5)") i,d_t_work(i)
2943       enddo
2944       Ek1=0.0d0
2945       ii=0
2946       do i=nnt,nct
2947 !        if (itype(i,1).eq.10) then
2948          mnum=molnum(i)
2949           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
2950         iti=iabs(itype(i,mnum))
2951 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2952          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2953           .or.(mnum.eq.5)) then
2954           j=ii+3
2955         else
2956           j=ii+6
2957         endif
2958         if (i.lt.nct) then
2959           do k=1,3
2960 !            write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2961             Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2962             0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2963           enddo
2964         endif
2965          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2966           .and.(mnum.ne.5)) ii=ii+3
2967         write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2968         write (iout,*) "ii",ii
2969         do k=1,3
2970           ii=ii+1
2971           write (iout,*) "k",k," ii",ii,"EK1",EK1
2972           if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2973           .and.(mnum.ne.5))&
2974            Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
2975           Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2976         enddo
2977         write (iout,*) "i",i," ii",ii
2978       enddo
2979       write (iout,*) "Ek from d_t_work",Ek1
2980       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2981 #endif
2982       if(allocated(DDU1)) deallocate(DDU1)
2983       if(allocated(DDU2)) deallocate(DDU2)
2984       if(allocated(DL2)) deallocate(DL2)
2985       if(allocated(DL1)) deallocate(DL1)
2986       if(allocated(xsolv)) deallocate(xsolv)
2987       if(allocated(DML)) deallocate(DML)
2988       if(allocated(rs)) deallocate(rs)
2989 #ifdef DEBUG
2990       if(allocated(matold)) deallocate(matold)
2991       if(allocated(rsold)) deallocate(rsold)
2992 #endif
2993       ind=1
2994       do j=nnt,nct
2995         do k=1,3
2996           d_t(k,j)=d_t_work(ind)
2997           ind=ind+1
2998         enddo
2999          mnum=molnum(j)
3000          if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3001           .and.(mnum.ne.5)) then
3002           do k=1,3
3003             d_t(k,j+nres)=d_t_work(ind)
3004             ind=ind+1
3005           enddo
3006         end if
3007       enddo
3008 #ifdef DEBUG
3009       write (iout,*) "Random velocities in the Calpha,SC space"
3010       do i=nnt,nct
3011         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3012       enddo
3013       do i=nnt,nct
3014         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3015       enddo
3016 #endif
3017       do j=1,3
3018         d_t(j,0)=d_t(j,nnt)
3019       enddo
3020       do i=nnt,nct
3021 !        if (itype(i,1).eq.10) then
3022          mnum=molnum(i)
3023          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3024           .or.(mnum.eq.5)) then
3025           do j=1,3
3026             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3027           enddo
3028         else
3029           do j=1,3
3030             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3031             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3032           enddo
3033         end if
3034       enddo
3035 #ifdef DEBUG
3036       write (iout,*)"Random velocities in the virtual-bond-vector space"
3037       do i=nnt,nct-1
3038         write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3039       enddo
3040       do i=nnt,nct
3041         write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3042       enddo
3043       call kinetic(Ek1)
3044       write (iout,*) "Ek from d_t_work",Ek1
3045       write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3046 #endif
3047 #else
3048       do k=0,2
3049         do i=1,dimen
3050           ind=(i-1)*3+k+1
3051           d_t_work(ind)=0.0d0
3052           do j=1,dimen
3053             d_t_work(ind)=d_t_work(ind) &
3054                             +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3055           enddo
3056 !          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3057 !          call flush(iout)
3058         enddo
3059       enddo
3060 ! Transfer to the d_t vector
3061       do j=1,3
3062         d_t(j,0)=d_t_work(j)
3063       enddo
3064       ind=3
3065       do i=nnt,nct-1
3066         do j=1,3
3067           ind=ind+1
3068           d_t(j,i)=d_t_work(ind)
3069         enddo
3070       enddo
3071       do i=nnt,nct
3072          mnum=molnum(i)
3073 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3074          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3075           .and.(mnum.ne.5)) then
3076           do j=1,3
3077             ind=ind+1
3078             d_t(j,i+nres)=d_t_work(ind)
3079           enddo
3080         endif
3081       enddo
3082 #endif
3083 !      call kinetic(EK)
3084 !      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3085 !        2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3086 !      call flush(iout)
3087 !      write(iout,*) "end init MD"
3088       return
3089       end subroutine random_vel
3090 !-----------------------------------------------------------------------------
3091 #ifndef LANG0
3092       subroutine sd_verlet_p_setup
3093 ! Sets up the parameters of stochastic Verlet algorithm
3094 !      implicit real*8 (a-h,o-z)
3095 !      include 'DIMENSIONS'
3096       use control, only: tcpu
3097       use control_data
3098 #ifdef MPI
3099       include 'mpif.h'
3100 #endif
3101 !      include 'COMMON.CONTROL'
3102 !      include 'COMMON.VAR'
3103 !      include 'COMMON.MD'
3104 !#ifndef LANG0
3105 !      include 'COMMON.LANGEVIN'
3106 !#else
3107 !      include 'COMMON.LANGEVIN.lang0'
3108 !#endif
3109 !      include 'COMMON.CHAIN'
3110 !      include 'COMMON.DERIV'
3111 !      include 'COMMON.GEO'
3112 !      include 'COMMON.LOCAL'
3113 !      include 'COMMON.INTERACT'
3114 !      include 'COMMON.IOUNITS'
3115 !      include 'COMMON.NAMES'
3116 !      include 'COMMON.TIME1'
3117       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
3118       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3119       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3120        prand_vec,vrand_vec1,vrand_vec2  !(MAXRES6) maxres6=6*maxres
3121       logical :: lprn = .false.
3122       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3123       real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3124                  gdt9,psig,tt0
3125       integer :: i,maxres2
3126 #ifdef MPI
3127       tt0 = MPI_Wtime()
3128 #else
3129       tt0 = tcpu()
3130 #endif
3131 !
3132 ! AL 8/17/04 Code adapted from tinker
3133 !
3134 ! Get the frictional and random terms for stochastic dynamics in the
3135 ! eigenspace of mass-scaled UNRES friction matrix
3136 !
3137       maxres2=2*nres
3138       do i = 1, dimen
3139             gdt = fricgam(i) * d_time
3140 !
3141 ! Stochastic dynamics reduces to simple MD for zero friction
3142 !
3143             if (gdt .le. zero) then
3144                pfric_vec(i) = 1.0d0
3145                vfric_vec(i) = d_time
3146                afric_vec(i) = 0.5d0 * d_time * d_time
3147                prand_vec(i) = 0.0d0
3148                vrand_vec1(i) = 0.0d0
3149                vrand_vec2(i) = 0.0d0
3150 !
3151 ! Analytical expressions when friction coefficient is large
3152 !
3153             else
3154                if (gdt .ge. gdt_radius) then
3155                   egdt = dexp(-gdt)
3156                   pfric_vec(i) = egdt
3157                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3158                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3159                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3160                   vterm = 1.0d0 - egdt**2
3161                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3162 !
3163 ! Use series expansions when friction coefficient is small
3164 !
3165                else
3166                   gdt2 = gdt * gdt
3167                   gdt3 = gdt * gdt2
3168                   gdt4 = gdt2 * gdt2
3169                   gdt5 = gdt2 * gdt3
3170                   gdt6 = gdt3 * gdt3
3171                   gdt7 = gdt3 * gdt4
3172                   gdt8 = gdt4 * gdt4
3173                   gdt9 = gdt4 * gdt5
3174                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3175                                 - gdt5/120.0d0 + gdt6/720.0d0 &
3176                                 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3177                                 - gdt9/362880.0d0) / fricgam(i)**2
3178                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3179                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3180                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3181                              + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3182                              + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3183                              + 127.0d0*gdt9/90720.0d0
3184                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3185                              - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3186                              - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3187                              - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3188                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3189                              - 17.0d0*gdt2/1280.0d0 &
3190                              + 17.0d0*gdt3/6144.0d0 &
3191                              + 40967.0d0*gdt4/34406400.0d0 &
3192                              - 57203.0d0*gdt5/275251200.0d0 &
3193                              - 1429487.0d0*gdt6/13212057600.0d0)
3194                end if
3195 !
3196 ! Compute the scaling factors of random terms for the nonzero friction case
3197 !
3198                ktm = 0.5d0*d_time/fricgam(i)
3199                psig = dsqrt(ktm*pterm) / fricgam(i)
3200                vsig = dsqrt(ktm*vterm)
3201                rhoc = dsqrt(1.0d0 - rho*rho)
3202                prand_vec(i) = psig
3203                vrand_vec1(i) = vsig * rho
3204                vrand_vec2(i) = vsig * rhoc
3205             end if
3206       end do
3207       if (lprn) then
3208       write (iout,*) &
3209         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3210         " vrand_vec2"
3211       do i=1,dimen
3212         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3213             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3214       enddo
3215       endif
3216 !
3217 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3218 !
3219 #ifndef   LANG0
3220       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3221       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3222       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3223       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3224       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3225       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3226 #endif
3227 #ifdef MPI
3228       t_sdsetup=t_sdsetup+MPI_Wtime()
3229 #else
3230       t_sdsetup=t_sdsetup+tcpu()-tt0
3231 #endif
3232       return
3233       end subroutine sd_verlet_p_setup
3234 !-----------------------------------------------------------------------------
3235       subroutine eigtransf1(n,ndim,ab,d,c)
3236
3237 !el      implicit none
3238       integer :: n,ndim
3239       real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3240       integer :: i,j,k
3241       do i=1,n
3242         do j=1,n
3243           c(i,j)=0.0d0
3244           do k=1,n
3245             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3246           enddo
3247         enddo
3248       enddo
3249       return
3250       end subroutine eigtransf1
3251 !-----------------------------------------------------------------------------
3252       subroutine eigtransf(n,ndim,a,b,d,c)
3253
3254 !el      implicit none
3255       integer :: n,ndim
3256       real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3257       integer :: i,j,k
3258       do i=1,n
3259         do j=1,n
3260           c(i,j)=0.0d0
3261           do k=1,n
3262             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3263           enddo
3264         enddo
3265       enddo
3266       return
3267       end subroutine eigtransf
3268 !-----------------------------------------------------------------------------
3269       subroutine sd_verlet1
3270
3271 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3272       use energy_data
3273 !      implicit real*8 (a-h,o-z)
3274 !      include 'DIMENSIONS'
3275 !      include 'COMMON.CONTROL'
3276 !      include 'COMMON.VAR'
3277 !      include 'COMMON.MD'
3278 !#ifndef LANG0
3279 !      include 'COMMON.LANGEVIN'
3280 !#else
3281 !      include 'COMMON.LANGEVIN.lang0'
3282 !#endif
3283 !      include 'COMMON.CHAIN'
3284 !      include 'COMMON.DERIV'
3285 !      include 'COMMON.GEO'
3286 !      include 'COMMON.LOCAL'
3287 !      include 'COMMON.INTERACT'
3288 !      include 'COMMON.IOUNITS'
3289 !      include 'COMMON.NAMES'
3290 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3291 !el      common /stochcalc/ stochforcvec
3292       logical :: lprn = .false.
3293       real(kind=8) :: ddt1,ddt2
3294       integer :: i,j,ind,inres
3295
3296 !      write (iout,*) "dc_old"
3297 !      do i=0,nres
3298 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3299 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3300 !      enddo
3301       do j=1,3
3302         dc_work(j)=dc_old(j,0)
3303         d_t_work(j)=d_t_old(j,0)
3304         d_a_work(j)=d_a_old(j,0)
3305       enddo
3306       ind=3
3307       do i=nnt,nct-1
3308         do j=1,3
3309           dc_work(ind+j)=dc_old(j,i)
3310           d_t_work(ind+j)=d_t_old(j,i)
3311           d_a_work(ind+j)=d_a_old(j,i)
3312         enddo
3313         ind=ind+3
3314       enddo
3315       do i=nnt,nct
3316          mnum=molnum(i)
3317 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3318          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3319           .and.(mnum.ne.5)) then
3320           do j=1,3
3321             dc_work(ind+j)=dc_old(j,i+nres)
3322             d_t_work(ind+j)=d_t_old(j,i+nres)
3323             d_a_work(ind+j)=d_a_old(j,i+nres)
3324           enddo
3325           ind=ind+3
3326         endif
3327       enddo
3328 #ifndef LANG0
3329       if (lprn) then
3330       write (iout,*) &
3331         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3332         " vrand_mat2"
3333       do i=1,dimen
3334         do j=1,dimen
3335           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3336             vfric_mat(i,j),afric_mat(i,j),&
3337             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3338         enddo
3339       enddo
3340       endif
3341       do i=1,dimen
3342         ddt1=0.0d0
3343         ddt2=0.0d0
3344         do j=1,dimen
3345           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3346             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3347           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3348           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3349         enddo
3350         d_t_work_new(i)=ddt1+0.5d0*ddt2
3351         d_t_work(i)=ddt1+ddt2
3352       enddo
3353 #endif
3354       do j=1,3
3355         dc(j,0)=dc_work(j)
3356         d_t(j,0)=d_t_work(j)
3357       enddo
3358       ind=3
3359       do i=nnt,nct-1
3360         do j=1,3
3361           dc(j,i)=dc_work(ind+j)
3362           d_t(j,i)=d_t_work(ind+j)
3363         enddo
3364         ind=ind+3
3365       enddo
3366       do i=nnt,nct
3367          mnum=molnum(i)
3368 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3369          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3370           .and.(mnum.ne.5)) then
3371           inres=i+nres
3372           do j=1,3
3373             dc(j,inres)=dc_work(ind+j)
3374             d_t(j,inres)=d_t_work(ind+j)
3375           enddo
3376           ind=ind+3
3377         endif
3378       enddo
3379       return
3380       end subroutine sd_verlet1
3381 !-----------------------------------------------------------------------------
3382       subroutine sd_verlet2
3383
3384 !  Calculating the adjusted velocities for accelerations
3385       use energy_data
3386 !      implicit real*8 (a-h,o-z)
3387 !      include 'DIMENSIONS'
3388 !      include 'COMMON.CONTROL'
3389 !      include 'COMMON.VAR'
3390 !      include 'COMMON.MD'
3391 !#ifndef LANG0
3392 !      include 'COMMON.LANGEVIN'
3393 !#else
3394 !      include 'COMMON.LANGEVIN.lang0'
3395 !#endif
3396 !      include 'COMMON.CHAIN'
3397 !      include 'COMMON.DERIV'
3398 !      include 'COMMON.GEO'
3399 !      include 'COMMON.LOCAL'
3400 !      include 'COMMON.INTERACT'
3401 !      include 'COMMON.IOUNITS'
3402 !      include 'COMMON.NAMES'
3403 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3404        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3405 !el      common /stochcalc/ stochforcvec
3406 !
3407       real(kind=8) :: ddt1,ddt2
3408       integer :: i,j,ind,inres
3409 ! Compute the stochastic forces which contribute to velocity change
3410 !
3411       call stochastic_force(stochforcvecV)
3412
3413 #ifndef LANG0
3414       do i=1,dimen
3415         ddt1=0.0d0
3416         ddt2=0.0d0
3417         do j=1,dimen
3418           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3419           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3420            vrand_mat2(i,j)*stochforcvecV(j)
3421         enddo
3422         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3423       enddo
3424 #endif
3425       do j=1,3
3426         d_t(j,0)=d_t_work(j)
3427       enddo
3428       ind=3
3429       do i=nnt,nct-1
3430         do j=1,3
3431           d_t(j,i)=d_t_work(ind+j)
3432         enddo
3433         ind=ind+3
3434       enddo
3435       do i=nnt,nct
3436          mnum=molnum(i)
3437 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3438          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3439           .and.(mnum.ne.5)) then
3440           inres=i+nres
3441           do j=1,3
3442             d_t(j,inres)=d_t_work(ind+j)
3443           enddo
3444           ind=ind+3
3445         endif
3446       enddo
3447       return
3448       end subroutine sd_verlet2
3449 !-----------------------------------------------------------------------------
3450       subroutine sd_verlet_ciccotti_setup
3451
3452 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3453 ! version
3454 !      implicit real*8 (a-h,o-z)
3455 !      include 'DIMENSIONS'
3456       use control, only: tcpu
3457       use control_data
3458 #ifdef MPI
3459       include 'mpif.h'
3460 #endif
3461 !      include 'COMMON.CONTROL'
3462 !      include 'COMMON.VAR'
3463 !      include 'COMMON.MD'
3464 !#ifndef LANG0
3465 !      include 'COMMON.LANGEVIN'
3466 !#else
3467 !      include 'COMMON.LANGEVIN.lang0'
3468 !#endif
3469 !      include 'COMMON.CHAIN'
3470 !      include 'COMMON.DERIV'
3471 !      include 'COMMON.GEO'
3472 !      include 'COMMON.LOCAL'
3473 !      include 'COMMON.INTERACT'
3474 !      include 'COMMON.IOUNITS'
3475 !      include 'COMMON.NAMES'
3476 !      include 'COMMON.TIME1'
3477       real(kind=8),dimension(6*nres) :: emgdt   !(MAXRES6) maxres6=6*maxres
3478       real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3479       real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3480         prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3481       logical :: lprn = .false.
3482       real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3483       real(kind=8) :: ktm,gdt,egdt,tt0
3484       integer :: i,maxres2
3485 #ifdef MPI
3486       tt0 = MPI_Wtime()
3487 #else
3488       tt0 = tcpu()
3489 #endif
3490 !
3491 ! AL 8/17/04 Code adapted from tinker
3492 !
3493 ! Get the frictional and random terms for stochastic dynamics in the
3494 ! eigenspace of mass-scaled UNRES friction matrix
3495 !
3496       maxres2=2*nres
3497       do i = 1, dimen
3498             write (iout,*) "i",i," fricgam",fricgam(i)
3499             gdt = fricgam(i) * d_time
3500 !
3501 ! Stochastic dynamics reduces to simple MD for zero friction
3502 !
3503             if (gdt .le. zero) then
3504                pfric_vec(i) = 1.0d0
3505                vfric_vec(i) = d_time
3506                afric_vec(i) = 0.5d0*d_time*d_time
3507                prand_vec(i) = afric_vec(i)
3508                vrand_vec2(i) = vfric_vec(i)
3509 !
3510 ! Analytical expressions when friction coefficient is large
3511 !
3512             else
3513                egdt = dexp(-gdt)
3514                pfric_vec(i) = egdt
3515                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3516                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3517                prand_vec(i) = afric_vec(i)
3518                vrand_vec2(i) = vfric_vec(i)
3519 !
3520 ! Compute the scaling factors of random terms for the nonzero friction case
3521 !
3522 !               ktm = 0.5d0*d_time/fricgam(i)
3523 !               psig = dsqrt(ktm*pterm) / fricgam(i)
3524 !               vsig = dsqrt(ktm*vterm)
3525 !               prand_vec(i) = psig*afric_vec(i)
3526 !               vrand_vec2(i) = vsig*vfric_vec(i)
3527             end if
3528       end do
3529       if (lprn) then
3530       write (iout,*) &
3531         "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3532         " vrand_vec2"
3533       do i=1,dimen
3534         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3535             afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3536       enddo
3537       endif
3538 !
3539 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3540 !
3541       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3542       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3543       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3544       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3545       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3546 #ifdef MPI
3547       t_sdsetup=t_sdsetup+MPI_Wtime()
3548 #else
3549       t_sdsetup=t_sdsetup+tcpu()-tt0
3550 #endif
3551       return
3552       end subroutine sd_verlet_ciccotti_setup
3553 !-----------------------------------------------------------------------------
3554       subroutine sd_verlet1_ciccotti
3555
3556 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3557 !      implicit real*8 (a-h,o-z)
3558       use energy_data
3559 !      include 'DIMENSIONS'
3560 #ifdef MPI
3561       include 'mpif.h'
3562 #endif
3563 !      include 'COMMON.CONTROL'
3564 !      include 'COMMON.VAR'
3565 !      include 'COMMON.MD'
3566 !#ifndef LANG0
3567 !      include 'COMMON.LANGEVIN'
3568 !#else
3569 !      include 'COMMON.LANGEVIN.lang0'
3570 !#endif
3571 !      include 'COMMON.CHAIN'
3572 !      include 'COMMON.DERIV'
3573 !      include 'COMMON.GEO'
3574 !      include 'COMMON.LOCAL'
3575 !      include 'COMMON.INTERACT'
3576 !      include 'COMMON.IOUNITS'
3577 !      include 'COMMON.NAMES'
3578 !el      real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3579 !el      common /stochcalc/ stochforcvec
3580       logical :: lprn = .false.
3581       real(kind=8) :: ddt1,ddt2
3582       integer :: i,j,ind,inres
3583 !      write (iout,*) "dc_old"
3584 !      do i=0,nres
3585 !        write (iout,'(i5,3f10.5,5x,3f10.5)')
3586 !     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3587 !      enddo
3588       do j=1,3
3589         dc_work(j)=dc_old(j,0)
3590         d_t_work(j)=d_t_old(j,0)
3591         d_a_work(j)=d_a_old(j,0)
3592       enddo
3593       ind=3
3594       do i=nnt,nct-1
3595         do j=1,3
3596           dc_work(ind+j)=dc_old(j,i)
3597           d_t_work(ind+j)=d_t_old(j,i)
3598           d_a_work(ind+j)=d_a_old(j,i)
3599         enddo
3600         ind=ind+3
3601       enddo
3602       do i=nnt,nct
3603         if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3604           do j=1,3
3605             dc_work(ind+j)=dc_old(j,i+nres)
3606             d_t_work(ind+j)=d_t_old(j,i+nres)
3607             d_a_work(ind+j)=d_a_old(j,i+nres)
3608           enddo
3609           ind=ind+3
3610         endif
3611       enddo
3612
3613 #ifndef LANG0
3614       if (lprn) then
3615       write (iout,*) &
3616         "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3617         " vrand_mat2"
3618       do i=1,dimen
3619         do j=1,dimen
3620           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3621             vfric_mat(i,j),afric_mat(i,j),&
3622             prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3623         enddo
3624       enddo
3625       endif
3626       do i=1,dimen
3627         ddt1=0.0d0
3628         ddt2=0.0d0
3629         do j=1,dimen
3630           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3631             +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3632           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3633           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3634         enddo
3635         d_t_work_new(i)=ddt1+0.5d0*ddt2
3636         d_t_work(i)=ddt1+ddt2
3637       enddo
3638 #endif
3639       do j=1,3
3640         dc(j,0)=dc_work(j)
3641         d_t(j,0)=d_t_work(j)
3642       enddo
3643       ind=3
3644       do i=nnt,nct-1
3645         do j=1,3
3646           dc(j,i)=dc_work(ind+j)
3647           d_t(j,i)=d_t_work(ind+j)
3648         enddo
3649         ind=ind+3
3650       enddo
3651       do i=nnt,nct
3652          mnum=molnum(i)
3653 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3654          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3655           .and.(mnum.ne.5)) then
3656           inres=i+nres
3657           do j=1,3
3658             dc(j,inres)=dc_work(ind+j)
3659             d_t(j,inres)=d_t_work(ind+j)
3660           enddo
3661           ind=ind+3
3662         endif
3663       enddo
3664       return
3665       end subroutine sd_verlet1_ciccotti
3666 !-----------------------------------------------------------------------------
3667       subroutine sd_verlet2_ciccotti
3668
3669 !  Calculating the adjusted velocities for accelerations
3670       use energy_data
3671 !      implicit real*8 (a-h,o-z)
3672 !      include 'DIMENSIONS'
3673 !      include 'COMMON.CONTROL'
3674 !      include 'COMMON.VAR'
3675 !      include 'COMMON.MD'
3676 !#ifndef LANG0
3677 !      include 'COMMON.LANGEVIN'
3678 !#else
3679 !      include 'COMMON.LANGEVIN.lang0'
3680 !#endif
3681 !      include 'COMMON.CHAIN'
3682 !      include 'COMMON.DERIV'
3683 !      include 'COMMON.GEO'
3684 !      include 'COMMON.LOCAL'
3685 !      include 'COMMON.INTERACT'
3686 !      include 'COMMON.IOUNITS'
3687 !      include 'COMMON.NAMES'
3688 !el      real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV   !(MAXRES6) maxres6=6*maxres
3689        real(kind=8),dimension(6*nres) :: stochforcvecV  !(MAXRES6) maxres6=6*maxres
3690 !el      common /stochcalc/ stochforcvec
3691       real(kind=8) :: ddt1,ddt2
3692       integer :: i,j,ind,inres
3693 !
3694 ! Compute the stochastic forces which contribute to velocity change
3695 !
3696       call stochastic_force(stochforcvecV)
3697 #ifndef LANG0
3698       do i=1,dimen
3699         ddt1=0.0d0
3700         ddt2=0.0d0
3701         do j=1,dimen
3702
3703           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3704 !          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3705           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3706         enddo
3707         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3708       enddo
3709 #endif
3710       do j=1,3
3711         d_t(j,0)=d_t_work(j)
3712       enddo
3713       ind=3
3714       do i=nnt,nct-1
3715         do j=1,3
3716           d_t(j,i)=d_t_work(ind+j)
3717         enddo
3718         ind=ind+3
3719       enddo
3720       do i=nnt,nct
3721          mnum=molnum(i)
3722          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3723           .and.(mnum.ne.5))
3724 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3725           inres=i+nres
3726           do j=1,3
3727             d_t(j,inres)=d_t_work(ind+j)
3728           enddo
3729           ind=ind+3
3730         endif
3731       enddo
3732       return
3733       end subroutine sd_verlet2_ciccotti
3734 #endif
3735 !-----------------------------------------------------------------------------
3736 ! moments.f
3737 !-----------------------------------------------------------------------------
3738       subroutine inertia_tensor
3739
3740 ! Calculating the intertia tensor for the entire protein in order to
3741 ! remove the perpendicular components of velocity matrix which cause
3742 ! the molecule to rotate.
3743       use comm_gucio
3744       use energy_data
3745 !       implicit real*8 (a-h,o-z)
3746 !       include 'DIMENSIONS'
3747 !       include 'COMMON.CONTROL'
3748 !       include 'COMMON.VAR'
3749 !       include 'COMMON.MD'
3750 !       include 'COMMON.CHAIN'
3751 !       include 'COMMON.DERIV'
3752 !       include 'COMMON.GEO'
3753 !       include 'COMMON.LOCAL'
3754 !       include 'COMMON.INTERACT'
3755 !       include 'COMMON.IOUNITS'
3756 !       include 'COMMON.NAMES'
3757
3758       real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3759       real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3760       real(kind=8) :: M_SC,mag,mag2,M_PEP
3761       real(kind=8),dimension(3,0:nres) :: vpp   !(3,0:MAXRES)
3762       real(kind=8),dimension(3) :: vs_p,pp,incr,v
3763       real(kind=8),dimension(3,3) :: pr1,pr2
3764
3765 !el      common /gucio/ cm
3766       integer :: iti,inres,i,j,k,mnum
3767         do i=1,3
3768            do j=1,3
3769              Im(i,j)=0.0d0
3770              pr1(i,j)=0.0d0
3771              pr2(i,j)=0.0d0
3772            enddo
3773            L(i)=0.0d0
3774            cm(i)=0.0d0
3775            vrot(i)=0.0d0
3776         enddo
3777 !   calculating the center of the mass of the protein
3778         M_PEP=0.0d0
3779         do i=nnt,nct-1
3780           mnum=molnum(i)
3781           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3782           if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3783           M_PEP=M_PEP+mp(mnum)
3784           do j=1,3
3785             cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3786           enddo
3787         enddo
3788 !        do j=1,3
3789 !         cm(j)=mp(1)*cm(j)
3790 !        enddo
3791         M_SC=0.0d0
3792         do i=nnt,nct
3793            mnum=molnum(i)
3794            if (mnum.eq.5) cycle
3795            iti=iabs(itype(i,mnum))
3796            M_SC=M_SC+msc(iabs(iti),mnum)
3797            inres=i+nres
3798            do j=1,3
3799             cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3800            enddo
3801         enddo
3802         do j=1,3
3803           cm(j)=cm(j)/(M_SC+M_PEP)
3804         enddo
3805
3806         do i=nnt,nct-1
3807            mnum=molnum(i)
3808           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3809           do j=1,3
3810             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3811           enddo
3812           Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3813           Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3814           Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3815           Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3816           Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3817           Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3818         enddo
3819
3820         do i=nnt,nct
3821            mnum=molnum(i)
3822            iti=iabs(itype(i,mnum))
3823            if (mnum.eq.5) cycle
3824            inres=i+nres
3825            do j=1,3
3826              pr(j)=c(j,inres)-cm(j)
3827            enddo
3828           Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3829           Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3830           Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3831           Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3832           Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3833           Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3834         enddo
3835
3836         do i=nnt,nct-1
3837            mnum=molnum(i)
3838           Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3839           vbld(i+1)*vbld(i+1)*0.25d0
3840           Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3841           vbld(i+1)*vbld(i+1)*0.25d0
3842           Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3843           vbld(i+1)*vbld(i+1)*0.25d0
3844           Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3845           vbld(i+1)*vbld(i+1)*0.25d0
3846           Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3847           vbld(i+1)*vbld(i+1)*0.25d0
3848           Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3849           vbld(i+1)*vbld(i+1)*0.25d0
3850         enddo
3851
3852
3853         do i=nnt,nct
3854               mnum=molnum(i)
3855 !         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3856          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3857           .and.(mnum.ne.5)) then
3858            iti=iabs(itype(i,mnum))
3859            inres=i+nres
3860           Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3861           dc_norm(1,inres))*vbld(inres)*vbld(inres)
3862           Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3863           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3864           Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3865           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3866           Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3867           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3868           Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3869           dc_norm(2,inres))*vbld(inres)*vbld(inres)
3870           Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3871           dc_norm(3,inres))*vbld(inres)*vbld(inres)
3872          endif
3873         enddo
3874
3875         call angmom(cm,L)
3876 !        write(iout,*) "The angular momentum before adjustment:"
3877 !        write(iout,*) (L(j),j=1,3)
3878
3879         Im(2,1)=Im(1,2)
3880         Im(3,1)=Im(1,3)
3881         Im(3,2)=Im(2,3)
3882
3883 !  Copying the Im matrix for the djacob subroutine
3884         do i=1,3
3885           do j=1,3
3886             Imcp(i,j)=Im(i,j)
3887             Id(i,j)=0.0d0
3888           enddo
3889         enddo
3890
3891 !   Finding the eigenvectors and eignvalues of the inertia tensor
3892        call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3893 !       write (iout,*) "Eigenvalues & Eigenvectors"
3894 !       write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3895 !       write (iout,*)
3896 !       do i=1,3
3897 !         write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3898 !       enddo
3899 !   Constructing the diagonalized matrix
3900        do i=1,3
3901          if (dabs(eigval(i)).gt.1.0d-15) then
3902            Id(i,i)=1.0d0/eigval(i)
3903          else
3904            Id(i,i)=0.0d0
3905          endif
3906        enddo
3907         do i=1,3
3908            do j=1,3
3909               Imcp(i,j)=eigvec(j,i)
3910            enddo
3911         enddo
3912         do i=1,3
3913            do j=1,3
3914               do k=1,3
3915                  pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3916               enddo
3917            enddo
3918         enddo
3919         do i=1,3
3920            do j=1,3
3921               do k=1,3
3922                  pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3923               enddo
3924            enddo
3925         enddo
3926 !  Calculating the total rotational velocity of the molecule
3927        do i=1,3
3928          do j=1,3
3929            vrot(i)=vrot(i)+pr2(i,j)*L(j)
3930          enddo
3931        enddo
3932 !   Resetting the velocities
3933        do i=nnt,nct-1
3934          call vecpr(vrot(1),dc(1,i),vp)
3935          do j=1,3
3936            d_t(j,i)=d_t(j,i)-vp(j)
3937           enddo
3938         enddo
3939         do i=nnt,nct
3940               mnum=molnum(i)
3941          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3942           .and.(mnum.ne.5)) then
3943 !         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3944 !       if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3945            inres=i+nres
3946            call vecpr(vrot(1),dc(1,inres),vp)
3947            do j=1,3
3948              d_t(j,inres)=d_t(j,inres)-vp(j)
3949            enddo
3950         endif
3951        enddo
3952        call angmom(cm,L)
3953 !       write(iout,*) "The angular momentum after adjustment:"
3954 !       write(iout,*) (L(j),j=1,3)
3955
3956       return
3957       end subroutine inertia_tensor
3958 !-----------------------------------------------------------------------------
3959       subroutine angmom(cm,L)
3960
3961       use energy_data
3962 !       implicit real*8 (a-h,o-z)
3963 !       include 'DIMENSIONS'
3964 !       include 'COMMON.CONTROL'
3965 !       include 'COMMON.VAR'
3966 !       include 'COMMON.MD'
3967 !       include 'COMMON.CHAIN'
3968 !       include 'COMMON.DERIV'
3969 !       include 'COMMON.GEO'
3970 !       include 'COMMON.LOCAL'
3971 !       include 'COMMON.INTERACT'
3972 !       include 'COMMON.IOUNITS'
3973 !       include 'COMMON.NAMES'
3974       real(kind=8) :: mscab
3975       real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3976       integer :: iti,inres,i,j,mnum
3977 !  Calculate the angular momentum
3978        do j=1,3
3979           L(j)=0.0d0
3980        enddo
3981        do j=1,3
3982           incr(j)=d_t(j,0)
3983        enddo
3984        do i=nnt,nct-1
3985           mnum=molnum(i)
3986           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3987           do j=1,3
3988             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3989           enddo
3990           do j=1,3
3991             v(j)=incr(j)+0.5d0*d_t(j,i)
3992           enddo
3993           do j=1,3
3994             incr(j)=incr(j)+d_t(j,i)
3995           enddo
3996           call vecpr(pr(1),v(1),vp)
3997           do j=1,3
3998             L(j)=L(j)+mp(mnum)*vp(j)
3999           enddo
4000           do j=1,3
4001              pr(j)=0.5d0*dc(j,i)
4002              pp(j)=0.5d0*d_t(j,i)
4003           enddo
4004          call vecpr(pr(1),pp(1),vp)
4005          do j=1,3
4006              L(j)=L(j)+Ip(mnum)*vp(j)
4007           enddo
4008         enddo
4009         do j=1,3
4010           incr(j)=d_t(j,0)
4011         enddo
4012         do i=nnt,nct
4013           mnum=molnum(i)
4014          iti=iabs(itype(i,mnum))
4015          inres=i+nres
4016         if (mnum.eq.5) then
4017          mscab=0.0d0
4018         else
4019          mscab=msc(iti,mnum)
4020         endif
4021          do j=1,3
4022            pr(j)=c(j,inres)-cm(j)
4023          enddo
4024          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4025           .and.(mnum.ne.5)) then
4026            do j=1,3
4027              v(j)=incr(j)+d_t(j,inres)
4028            enddo
4029          else
4030            do j=1,3
4031              v(j)=incr(j)
4032            enddo
4033          endif
4034          call vecpr(pr(1),v(1),vp)
4035 !         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4036 !           " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4037          do j=1,3
4038             L(j)=L(j)+mscab*vp(j)
4039          enddo
4040 !         write (iout,*) "L",(l(j),j=1,3)
4041          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4042           .and.(mnum.ne.5)) then
4043            do j=1,3
4044             v(j)=incr(j)+d_t(j,inres)
4045            enddo
4046            call vecpr(dc(1,inres),d_t(1,inres),vp)
4047            do j=1,3
4048              L(j)=L(j)+Isc(iti,mnum)*vp(j)
4049           enddo
4050          endif
4051          do j=1,3
4052              incr(j)=incr(j)+d_t(j,i)
4053          enddo
4054        enddo
4055       return
4056       end subroutine angmom
4057 !-----------------------------------------------------------------------------
4058       subroutine vcm_vel(vcm)
4059
4060       use energy_data
4061 !       implicit real*8 (a-h,o-z)
4062 !       include 'DIMENSIONS'
4063 !       include 'COMMON.VAR'
4064 !       include 'COMMON.MD'
4065 !       include 'COMMON.CHAIN'
4066 !       include 'COMMON.DERIV'
4067 !       include 'COMMON.GEO'
4068 !       include 'COMMON.LOCAL'
4069 !       include 'COMMON.INTERACT'
4070 !       include 'COMMON.IOUNITS'
4071        real(kind=8),dimension(3) :: vcm,vv
4072        real(kind=8) :: summas,amas
4073        integer :: i,j,mnum
4074
4075        do j=1,3
4076          vcm(j)=0.0d0
4077          vv(j)=d_t(j,0)
4078        enddo
4079        summas=0.0d0
4080        do i=nnt,nct
4081          mnum=molnum(i)
4082          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4083          if (i.lt.nct) then
4084            summas=summas+mp(mnum)
4085            do j=1,3
4086              vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4087            enddo
4088          endif
4089          if (mnum.ne.5) then
4090          amas=msc(iabs(itype(i,mnum)),mnum)
4091          else
4092          amas=0.0d0
4093          endif
4094          summas=summas+amas
4095          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4096           .and.(mnum.ne.5)) then
4097            do j=1,3
4098              vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4099            enddo
4100          else
4101            do j=1,3
4102              vcm(j)=vcm(j)+amas*vv(j)
4103            enddo
4104          endif
4105          do j=1,3
4106            vv(j)=vv(j)+d_t(j,i)
4107          enddo
4108        enddo
4109 !       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4110        do j=1,3
4111          vcm(j)=vcm(j)/summas
4112        enddo
4113       return
4114       end subroutine vcm_vel
4115 !-----------------------------------------------------------------------------
4116 ! rattle.F
4117 !-----------------------------------------------------------------------------
4118       subroutine rattle1
4119 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4120 ! AL 9/24/04
4121       use comm_przech
4122       use energy_data
4123 !      implicit real*8 (a-h,o-z)
4124 !      include 'DIMENSIONS'
4125 #ifdef RATTLE
4126 !      include 'COMMON.CONTROL'
4127 !      include 'COMMON.VAR'
4128 !      include 'COMMON.MD'
4129 !#ifndef LANG0
4130 !      include 'COMMON.LANGEVIN'
4131 !#else
4132 !      include 'COMMON.LANGEVIN.lang0'
4133 !#endif
4134 !      include 'COMMON.CHAIN'
4135 !      include 'COMMON.DERIV'
4136 !      include 'COMMON.GEO'
4137 !      include 'COMMON.LOCAL'
4138 !      include 'COMMON.INTERACT'
4139 !      include 'COMMON.IOUNITS'
4140 !      include 'COMMON.NAMES'
4141 !      include 'COMMON.TIME1'
4142 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4143 !el       gdc(3,2*nres,2*nres)
4144       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4145 !el      real(kind=8) :: Cmat(2*nres,2*nres)
4146       real(kind=8) :: x(2*nres),xcorr(3,2*nres)         !maxres2=2*maxres
4147 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4148 !el      common /przechowalnia/ nbond
4149       integer :: max_rattle = 5
4150       logical :: lprn = .false., lprn1 = .false., not_done
4151       real(kind=8) :: tol_rattle = 1.0d-5
4152
4153       integer :: ii,i,j,jj,l,ind,ind1,nres2
4154       nres2=2*nres
4155
4156 !el /common/ przechowalnia
4157
4158       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4159       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4160       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4161 !el--------
4162       if (lprn) write (iout,*) "RATTLE1"
4163       nbond=nct-nnt
4164       do i=nnt,nct
4165        mnum=molnum(i)
4166          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4167           .and.(mnum.ne.5)) nbond=nbond+1
4168       enddo
4169 ! Make a folded form of the Ginv-matrix
4170       ind=0
4171       ii=0
4172       do i=nnt,nct-1
4173         ii=ii+1
4174         do j=1,3
4175           ind=ind+1
4176           ind1=0
4177           jj=0
4178           do k=nnt,nct-1
4179             jj=jj+1
4180             do l=1,3
4181               ind1=ind1+1
4182               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4183             enddo
4184           enddo
4185           do k=nnt,nct
4186           mnum=molnum(k)
4187          if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4188           .and.(mnum.ne.5)) then
4189               jj=jj+1
4190               do l=1,3
4191                 ind1=ind1+1
4192                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4193               enddo
4194             endif
4195           enddo
4196         enddo
4197       enddo
4198       do i=nnt,nct
4199          mnum=molnum(i)
4200          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4201           .and.(mnum.ne.5))
4202           ii=ii+1
4203           do j=1,3
4204             ind=ind+1
4205             ind1=0
4206             jj=0
4207             do k=nnt,nct-1
4208               jj=jj+1
4209               do l=1,3
4210                 ind1=ind1+1
4211                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4212               enddo
4213             enddo
4214             do k=nnt,nct
4215               if (itype(k,1).ne.10) then
4216                 jj=jj+1
4217                 do l=1,3
4218                   ind1=ind1+1
4219                   if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4220                 enddo
4221               endif
4222             enddo
4223           enddo
4224         endif
4225       enddo
4226       if (lprn1) then
4227         write (iout,*) "Matrix GGinv"
4228         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4229       endif
4230       not_done=.true.
4231       iter=0
4232       do while (not_done)
4233       iter=iter+1
4234       if (iter.gt.max_rattle) then
4235         write (iout,*) "Error - too many iterations in RATTLE."
4236         stop
4237       endif
4238 ! Calculate the matrix C = GG**(-1) dC_old o dC
4239       ind1=0
4240       do i=nnt,nct-1
4241         ind1=ind1+1
4242         do j=1,3
4243           dC_uncor(j,ind1)=dC(j,i)
4244         enddo
4245       enddo
4246       do i=nnt,nct
4247         if (itype(i,1).ne.10) then
4248           ind1=ind1+1
4249           do j=1,3
4250             dC_uncor(j,ind1)=dC(j,i+nres)
4251           enddo
4252         endif
4253       enddo
4254       do i=1,nbond
4255         ind=0
4256         do k=nnt,nct-1
4257           ind=ind+1
4258           do j=1,3
4259             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4260           enddo
4261         enddo
4262         do k=nnt,nct
4263           if (itype(k,1).ne.10) then
4264             ind=ind+1
4265             do j=1,3
4266               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4267             enddo
4268           endif
4269         enddo
4270       enddo
4271 ! Calculate deviations from standard virtual-bond lengths
4272       ind=0
4273       do i=nnt,nct-1
4274         ind=ind+1
4275         x(ind)=vbld(i+1)**2-vbl**2
4276       enddo
4277       do i=nnt,nct
4278         if (itype(i,1).ne.10) then
4279           ind=ind+1
4280           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4281         endif
4282       enddo
4283       if (lprn) then
4284         write (iout,*) "Coordinates and violations"
4285         do i=1,nbond
4286           write(iout,'(i5,3f10.5,5x,e15.5)') &
4287            i,(dC_uncor(j,i),j=1,3),x(i)
4288         enddo
4289         write (iout,*) "Velocities and violations"
4290         ind=0
4291         do i=nnt,nct-1
4292           ind=ind+1
4293           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4294            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4295         enddo
4296         do i=nnt,nct
4297           mnum=molnum(i)
4298          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4299           .and.(mnum.ne.5)) then
4300
4301             ind=ind+1
4302             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4303              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4304              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4305           endif
4306         enddo
4307 !        write (iout,*) "gdc"
4308 !        do i=1,nbond
4309 !          write (iout,*) "i",i
4310 !          do j=1,nbond
4311 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4312 !          enddo
4313 !        enddo
4314       endif
4315       xmax=dabs(x(1))
4316       do i=2,nbond
4317         if (dabs(x(i)).gt.xmax) then
4318           xmax=dabs(x(i))
4319         endif
4320       enddo
4321       if (xmax.lt.tol_rattle) then
4322         not_done=.false.
4323         goto 100
4324       endif
4325 ! Calculate the matrix of the system of equations
4326       do i=1,nbond
4327         do j=1,nbond
4328           Cmat(i,j)=0.0d0
4329           do k=1,3
4330             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4331           enddo
4332         enddo
4333       enddo
4334       if (lprn1) then
4335         write (iout,*) "Matrix Cmat"
4336         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4337       endif
4338       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4339 ! Add constraint term to positions
4340       ind=0
4341       do i=nnt,nct-1
4342         ind=ind+1
4343         do j=1,3
4344           xx=0.0d0
4345           do ii=1,nbond
4346             xx = xx+x(ii)*gdc(j,ind,ii)
4347           enddo
4348           xx=0.5d0*xx
4349           dC(j,i)=dC(j,i)-xx
4350           d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4351         enddo
4352       enddo
4353       do i=nnt,nct
4354         if (itype(i,1).ne.10) then
4355           ind=ind+1
4356           do j=1,3
4357             xx=0.0d0
4358             do ii=1,nbond
4359               xx = xx+x(ii)*gdc(j,ind,ii)
4360             enddo
4361             xx=0.5d0*xx
4362             dC(j,i+nres)=dC(j,i+nres)-xx
4363             d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4364           enddo
4365         endif
4366       enddo
4367 ! Rebuild the chain using the new coordinates
4368       call chainbuild_cart
4369       if (lprn) then
4370         write (iout,*) "New coordinates, Lagrange multipliers,",&
4371         " and differences between actual and standard bond lengths"
4372         ind=0
4373         do i=nnt,nct-1
4374           ind=ind+1
4375           xx=vbld(i+1)**2-vbl**2
4376           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4377               i,(dC(j,i),j=1,3),x(ind),xx
4378         enddo
4379         do i=nnt,nct
4380          mnum=molnum(i)
4381          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4382           .and.(mnum.ne.5))
4383             ind=ind+1
4384             xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4385             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4386              i,(dC(j,i+nres),j=1,3),x(ind),xx
4387           endif
4388         enddo
4389         write (iout,*) "Velocities and violations"
4390         ind=0
4391         do i=nnt,nct-1
4392           ind=ind+1
4393           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4394            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4395         enddo
4396         do i=nnt,nct
4397           if (itype(i,1).ne.10) then
4398             ind=ind+1
4399             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4400              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4401              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4402           endif
4403         enddo
4404       endif
4405       enddo
4406   100 continue
4407       return
4408    10 write (iout,*) "Error - singularity in solving the system",&
4409        " of equations for Lagrange multipliers."
4410       stop
4411 #else
4412       write (iout,*) &
4413        "RATTLE inactive; use -DRATTLE switch at compile time."
4414       stop
4415 #endif
4416       end subroutine rattle1
4417 !-----------------------------------------------------------------------------
4418       subroutine rattle2
4419 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4420 ! AL 9/24/04
4421       use comm_przech
4422       use energy_data
4423 !      implicit real*8 (a-h,o-z)
4424 !      include 'DIMENSIONS'
4425 #ifdef RATTLE
4426 !      include 'COMMON.CONTROL'
4427 !      include 'COMMON.VAR'
4428 !      include 'COMMON.MD'
4429 !#ifndef LANG0
4430 !      include 'COMMON.LANGEVIN'
4431 !#else
4432 !      include 'COMMON.LANGEVIN.lang0'
4433 !#endif
4434 !      include 'COMMON.CHAIN'
4435 !      include 'COMMON.DERIV'
4436 !      include 'COMMON.GEO'
4437 !      include 'COMMON.LOCAL'
4438 !      include 'COMMON.INTERACT'
4439 !      include 'COMMON.IOUNITS'
4440 !      include 'COMMON.NAMES'
4441 !      include 'COMMON.TIME1'
4442 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4443 !el       gdc(3,2*nres,2*nres)
4444       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4445 !el       Cmat(2*nres,2*nres)
4446       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4447 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4448 !el      common /przechowalnia/ nbond
4449       integer :: max_rattle = 5
4450       logical :: lprn = .false., lprn1 = .false., not_done
4451       real(kind=8) :: tol_rattle = 1.0d-5
4452       integer :: nres2
4453       nres2=2*nres
4454
4455 !el /common/ przechowalnia
4456       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4457       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4458       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4459 !el--------
4460       if (lprn) write (iout,*) "RATTLE2"
4461       if (lprn) write (iout,*) "Velocity correction"
4462 ! Calculate the matrix G dC
4463       do i=1,nbond
4464         ind=0
4465         do k=nnt,nct-1
4466           ind=ind+1
4467           do j=1,3
4468             gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4469           enddo
4470         enddo
4471         do k=nnt,nct
4472          mnum=molnum(i)
4473          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4474           .and.(mnum.ne.5)) then
4475             ind=ind+1
4476             do j=1,3
4477               gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4478             enddo
4479           endif
4480         enddo
4481       enddo
4482 !      if (lprn) then
4483 !        write (iout,*) "gdc"
4484 !        do i=1,nbond
4485 !          write (iout,*) "i",i
4486 !          do j=1,nbond
4487 !            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4488 !          enddo
4489 !        enddo
4490 !      endif
4491 ! Calculate the matrix of the system of equations
4492       ind=0
4493       do i=nnt,nct-1
4494         ind=ind+1
4495         do j=1,nbond
4496           Cmat(ind,j)=0.0d0
4497           do k=1,3
4498             Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4499           enddo
4500         enddo
4501       enddo
4502       do i=nnt,nct
4503          mnum=molnum(i)
4504          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4505           .and.(mnum.ne.5)) then
4506           ind=ind+1
4507           do j=1,nbond
4508             Cmat(ind,j)=0.0d0
4509             do k=1,3
4510               Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4511             enddo
4512           enddo
4513         endif
4514       enddo
4515 ! Calculate the scalar product dC o d_t_new
4516       ind=0
4517       do i=nnt,nct-1
4518         ind=ind+1
4519         x(ind)=scalar(d_t(1,i),dC(1,i))
4520       enddo
4521       do i=nnt,nct
4522          mnum=molnum(i)
4523          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4524           .and.(mnum.ne.5)) then
4525           ind=ind+1
4526           x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4527         endif
4528       enddo
4529       if (lprn) then
4530         write (iout,*) "Velocities and violations"
4531         ind=0
4532         do i=nnt,nct-1
4533           ind=ind+1
4534           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4535            i,ind,(d_t(j,i),j=1,3),x(ind)
4536         enddo
4537         do i=nnt,nct
4538          mnum=molnum(i)
4539          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4540           .and.(mnum.ne.5)) then
4541             ind=ind+1
4542             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4543              i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4544           endif
4545         enddo
4546       endif
4547       xmax=dabs(x(1))
4548       do i=2,nbond
4549         if (dabs(x(i)).gt.xmax) then
4550           xmax=dabs(x(i))
4551         endif
4552       enddo
4553       if (xmax.lt.tol_rattle) then
4554         not_done=.false.
4555         goto 100
4556       endif
4557       if (lprn1) then
4558         write (iout,*) "Matrix Cmat"
4559         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4560       endif
4561       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4562 ! Add constraint term to velocities
4563       ind=0
4564       do i=nnt,nct-1
4565         ind=ind+1
4566         do j=1,3
4567           xx=0.0d0
4568           do ii=1,nbond
4569             xx = xx+x(ii)*gdc(j,ind,ii)
4570           enddo
4571           d_t(j,i)=d_t(j,i)-xx
4572         enddo
4573       enddo
4574       do i=nnt,nct
4575          mnum=molnum(i)
4576          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4577           .and.(mnum.ne.5)) then
4578           ind=ind+1
4579           do j=1,3
4580             xx=0.0d0
4581             do ii=1,nbond
4582               xx = xx+x(ii)*gdc(j,ind,ii)
4583             enddo
4584             d_t(j,i+nres)=d_t(j,i+nres)-xx
4585           enddo
4586         endif
4587       enddo
4588       if (lprn) then
4589         write (iout,*) &
4590           "New velocities, Lagrange multipliers violations"
4591         ind=0
4592         do i=nnt,nct-1
4593           ind=ind+1
4594           if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4595              i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4596         enddo
4597         do i=nnt,nct
4598          mnum=molnum(i)
4599          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4600           .and.(mnum.ne.5))
4601             ind=ind+1
4602             write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4603               i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4604               scalar(d_t(1,i+nres),dC(1,i+nres))
4605           endif
4606         enddo
4607       endif
4608   100 continue
4609       return
4610    10 write (iout,*) "Error - singularity in solving the system",&
4611        " of equations for Lagrange multipliers."
4612       stop
4613 #else
4614       write (iout,*) &
4615        "RATTLE inactive; use -DRATTLE option at compile time."
4616       stop
4617 #endif
4618       end subroutine rattle2
4619 !-----------------------------------------------------------------------------
4620       subroutine rattle_brown
4621 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4622 ! AL 9/24/04
4623       use comm_przech
4624       use energy_data
4625 !      implicit real*8 (a-h,o-z)
4626 !      include 'DIMENSIONS'
4627 #ifdef RATTLE
4628 !      include 'COMMON.CONTROL'
4629 !      include 'COMMON.VAR'
4630 !      include 'COMMON.MD'
4631 !#ifndef LANG0
4632 !      include 'COMMON.LANGEVIN'
4633 !#else
4634 !      include 'COMMON.LANGEVIN.lang0'
4635 !#endif
4636 !      include 'COMMON.CHAIN'
4637 !      include 'COMMON.DERIV'
4638 !      include 'COMMON.GEO'
4639 !      include 'COMMON.LOCAL'
4640 !      include 'COMMON.INTERACT'
4641 !      include 'COMMON.IOUNITS'
4642 !      include 'COMMON.NAMES'
4643 !      include 'COMMON.TIME1'
4644 !el      real(kind=8) :: gginv(2*nres,2*nres),&
4645 !el       gdc(3,2*nres,2*nres)
4646       real(kind=8) :: dC_uncor(3,2*nres)        !,&
4647 !el      real(kind=8) :: Cmat(2*nres,2*nres)
4648       real(kind=8) :: x(2*nres)         !maxres2=2*maxres
4649 !el      common /przechowalnia/ GGinv,gdc,Cmat,nbond
4650 !el      common /przechowalnia/ nbond
4651       integer :: max_rattle = 5
4652       logical :: lprn = .false., lprn1 = .false., not_done
4653       real(kind=8) :: tol_rattle = 1.0d-5
4654       integer :: nres2
4655       nres2=2*nres
4656
4657 !el /common/ przechowalnia
4658       if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4659       if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4660       if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4661 !el--------
4662
4663       if (lprn) write (iout,*) "RATTLE_BROWN"
4664       nbond=nct-nnt
4665       do i=nnt,nct
4666         if (itype(i,1).ne.10) nbond=nbond+1
4667       enddo
4668 ! Make a folded form of the Ginv-matrix
4669       ind=0
4670       ii=0
4671       do i=nnt,nct-1
4672         ii=ii+1
4673         do j=1,3
4674           ind=ind+1
4675           ind1=0
4676           jj=0
4677           do k=nnt,nct-1
4678             jj=jj+1
4679             do l=1,3
4680               ind1=ind1+1
4681               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4682             enddo
4683           enddo
4684           do k=nnt,nct
4685             if (itype(k,1).ne.10) then
4686               jj=jj+1
4687               do l=1,3
4688                 ind1=ind1+1
4689                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4690               enddo
4691             endif
4692           enddo
4693         enddo
4694       enddo
4695       do i=nnt,nct
4696         if (itype(i,1).ne.10) then
4697           ii=ii+1
4698           do j=1,3
4699             ind=ind+1
4700             ind1=0
4701             jj=0
4702             do k=nnt,nct-1
4703               jj=jj+1
4704               do l=1,3
4705                 ind1=ind1+1
4706                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4707               enddo
4708             enddo
4709             do k=nnt,nct
4710               if (itype(k,1).ne.10) then
4711                 jj=jj+1
4712                 do l=1,3
4713                   ind1=ind1+1
4714                   if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4715                 enddo
4716               endif
4717             enddo
4718           enddo
4719         endif
4720       enddo
4721       if (lprn1) then
4722         write (iout,*) "Matrix GGinv"
4723         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4724       endif
4725       not_done=.true.
4726       iter=0
4727       do while (not_done)
4728       iter=iter+1
4729       if (iter.gt.max_rattle) then
4730         write (iout,*) "Error - too many iterations in RATTLE."
4731         stop
4732       endif
4733 ! Calculate the matrix C = GG**(-1) dC_old o dC
4734       ind1=0
4735       do i=nnt,nct-1
4736         ind1=ind1+1
4737         do j=1,3
4738           dC_uncor(j,ind1)=dC(j,i)
4739         enddo
4740       enddo
4741       do i=nnt,nct
4742         if (itype(i,1).ne.10) then
4743           ind1=ind1+1
4744           do j=1,3
4745             dC_uncor(j,ind1)=dC(j,i+nres)
4746           enddo
4747         endif
4748       enddo
4749       do i=1,nbond
4750         ind=0
4751         do k=nnt,nct-1
4752           ind=ind+1
4753           do j=1,3
4754             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4755           enddo
4756         enddo
4757         do k=nnt,nct
4758           if (itype(k,1).ne.10) then
4759             ind=ind+1
4760             do j=1,3
4761               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4762             enddo
4763           endif
4764         enddo
4765       enddo
4766 ! Calculate deviations from standard virtual-bond lengths
4767       ind=0
4768       do i=nnt,nct-1
4769         ind=ind+1
4770         x(ind)=vbld(i+1)**2-vbl**2
4771       enddo
4772       do i=nnt,nct
4773         if (itype(i,1).ne.10) then
4774           ind=ind+1
4775           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4776         endif
4777       enddo
4778       if (lprn) then
4779         write (iout,*) "Coordinates and violations"
4780         do i=1,nbond
4781           write(iout,'(i5,3f10.5,5x,e15.5)') &
4782            i,(dC_uncor(j,i),j=1,3),x(i)
4783         enddo
4784         write (iout,*) "Velocities and violations"
4785         ind=0
4786         do i=nnt,nct-1
4787           ind=ind+1
4788           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4789            i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4790         enddo
4791         do i=nnt,nct
4792           if (itype(i,1).ne.10) then
4793             ind=ind+1
4794             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4795              i+nres,ind,(d_t(j,i+nres),j=1,3),&
4796              scalar(d_t(1,i+nres),dC_old(1,i+nres))
4797           endif
4798         enddo
4799         write (iout,*) "gdc"
4800         do i=1,nbond
4801           write (iout,*) "i",i
4802           do j=1,nbond
4803             write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4804           enddo
4805         enddo
4806       endif
4807       xmax=dabs(x(1))
4808       do i=2,nbond
4809         if (dabs(x(i)).gt.xmax) then
4810           xmax=dabs(x(i))
4811         endif
4812       enddo
4813       if (xmax.lt.tol_rattle) then
4814         not_done=.false.
4815         goto 100
4816       endif
4817 ! Calculate the matrix of the system of equations
4818       do i=1,nbond
4819         do j=1,nbond
4820           Cmat(i,j)=0.0d0
4821           do k=1,3
4822             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4823           enddo
4824         enddo
4825       enddo
4826       if (lprn1) then
4827         write (iout,*) "Matrix Cmat"
4828         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4829       endif
4830       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4831 ! Add constraint term to positions
4832       ind=0
4833       do i=nnt,nct-1
4834         ind=ind+1
4835         do j=1,3
4836           xx=0.0d0
4837           do ii=1,nbond
4838             xx = xx+x(ii)*gdc(j,ind,ii)
4839           enddo
4840           xx=-0.5d0*xx
4841           d_t(j,i)=d_t(j,i)+xx/d_time
4842           dC(j,i)=dC(j,i)+xx
4843         enddo
4844       enddo
4845       do i=nnt,nct
4846         if (itype(i,1).ne.10) then
4847           ind=ind+1
4848           do j=1,3
4849             xx=0.0d0
4850             do ii=1,nbond
4851               xx = xx+x(ii)*gdc(j,ind,ii)
4852             enddo
4853             xx=-0.5d0*xx
4854             d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4855             dC(j,i+nres)=dC(j,i+nres)+xx
4856           enddo
4857         endif
4858       enddo
4859 ! Rebuild the chain using the new coordinates
4860       call chainbuild_cart
4861       if (lprn) then
4862         write (iout,*) "New coordinates, Lagrange multipliers,",&
4863         " and differences between actual and standard bond lengths"
4864         ind=0
4865         do i=nnt,nct-1
4866           ind=ind+1
4867           xx=vbld(i+1)**2-vbl**2
4868           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4869               i,(dC(j,i),j=1,3),x(ind),xx
4870         enddo
4871         do i=nnt,nct
4872           if (itype(i,1).ne.10) then
4873             ind=ind+1
4874             xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4875             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4876              i,(dC(j,i+nres),j=1,3),x(ind),xx
4877           endif
4878         enddo
4879         write (iout,*) "Velocities and violations"
4880         ind=0
4881         do i=nnt,nct-1
4882           ind=ind+1
4883           write (iout,'(2i5,3f10.5,5x,e15.5)') &
4884            i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4885         enddo
4886         do i=nnt,nct
4887           if (itype(i,1).ne.10) then
4888             ind=ind+1
4889             write (iout,'(2i5,3f10.5,5x,e15.5)') &
4890              i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4891              scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4892           endif
4893         enddo
4894       endif
4895       enddo
4896   100 continue
4897       return
4898    10 write (iout,*) "Error - singularity in solving the system",&
4899        " of equations for Lagrange multipliers."
4900       stop
4901 #else
4902       write (iout,*) &
4903        "RATTLE inactive; use -DRATTLE option at compile time"
4904       stop
4905 #endif
4906       end subroutine rattle_brown
4907 !-----------------------------------------------------------------------------
4908 ! stochfric.F
4909 !-----------------------------------------------------------------------------
4910       subroutine friction_force
4911
4912       use energy_data
4913       use REMD_data
4914       use comm_syfek
4915 !      implicit real*8 (a-h,o-z)
4916 !      include 'DIMENSIONS'
4917 !      include 'COMMON.VAR'
4918 !      include 'COMMON.CHAIN'
4919 !      include 'COMMON.DERIV'
4920 !      include 'COMMON.GEO'
4921 !      include 'COMMON.LOCAL'
4922 !      include 'COMMON.INTERACT'
4923 !      include 'COMMON.MD'
4924 !#ifndef LANG0
4925 !      include 'COMMON.LANGEVIN'
4926 !#else
4927 !      include 'COMMON.LANGEVIN.lang0'
4928 !#endif
4929 !      include 'COMMON.IOUNITS'
4930 !el      real(kind=8),dimension(6*nres) :: gamvec       !(MAXRES6) maxres6=6*maxres
4931 !el      common /syfek/ gamvec
4932       real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4933 !el       ginvfric(2*nres,2*nres)       !maxres2=2*maxres
4934 !el      common /przechowalnia/ ginvfric
4935
4936       logical :: lprn = .false., checkmode = .false.
4937       integer :: i,j,ind,k,nres2,nres6,mnum
4938       nres2=2*nres
4939       nres6=6*nres
4940
4941       if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4942       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4943       do i=0,nres2
4944         do j=1,3
4945           friction(j,i)=0.0d0
4946         enddo
4947       enddo
4948
4949       do j=1,3
4950         d_t_work(j)=d_t(j,0)
4951       enddo
4952       ind=3
4953       do i=nnt,nct-1
4954         do j=1,3
4955           d_t_work(ind+j)=d_t(j,i)
4956         enddo
4957         ind=ind+3
4958       enddo
4959       do i=nnt,nct
4960         mnum=molnum(i)
4961         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4962         .and.(mnum.ne.5)) then
4963           do j=1,3
4964             d_t_work(ind+j)=d_t(j,i+nres)
4965           enddo
4966           ind=ind+3
4967         endif
4968       enddo
4969
4970       call fricmat_mult(d_t_work,fric_work)
4971
4972       if (.not.checkmode) return
4973
4974       if (lprn) then
4975         write (iout,*) "d_t_work and fric_work"
4976         do i=1,3*dimen
4977           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4978         enddo
4979       endif
4980       do j=1,3
4981         friction(j,0)=fric_work(j)
4982       enddo
4983       ind=3
4984       do i=nnt,nct-1
4985         do j=1,3
4986           friction(j,i)=fric_work(ind+j)
4987         enddo
4988         ind=ind+3
4989       enddo
4990       do i=nnt,nct
4991         mnum=molnum(i)
4992         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4993         .and.(mnum.ne.5)) then
4994 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4995           do j=1,3
4996             friction(j,i+nres)=fric_work(ind+j)
4997           enddo
4998           ind=ind+3
4999         endif
5000       enddo
5001       if (lprn) then
5002         write(iout,*) "Friction backbone"
5003         do i=0,nct-1
5004           write(iout,'(i5,3e15.5,5x,3e15.5)') &
5005            i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5006         enddo
5007         write(iout,*) "Friction side chain"
5008         do i=nnt,nct
5009           write(iout,'(i5,3e15.5,5x,3e15.5)') &
5010            i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5011         enddo
5012       endif
5013       if (lprn) then
5014         do j=1,3
5015           vv(j)=d_t(j,0)
5016         enddo
5017         do i=nnt,nct
5018           do j=1,3
5019             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5020             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5021             vv(j)=vv(j)+d_t(j,i)
5022           enddo
5023         enddo
5024         write (iout,*) "vvtot backbone and sidechain"
5025         do i=nnt,nct
5026           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5027            (vvtot(j,i+nres),j=1,3)
5028         enddo
5029         ind=0
5030         do i=nnt,nct-1
5031           do j=1,3
5032             v_work(ind+j)=vvtot(j,i)
5033           enddo
5034           ind=ind+3
5035         enddo
5036         do i=nnt,nct
5037           do j=1,3
5038             v_work(ind+j)=vvtot(j,i+nres)
5039           enddo
5040           ind=ind+3
5041         enddo
5042         write (iout,*) "v_work gamvec and site-based friction forces"
5043         do i=1,dimen1
5044           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5045             gamvec(i)*v_work(i)
5046         enddo
5047 !        do i=1,dimen
5048 !          fric_work1(i)=0.0d0
5049 !          do j=1,dimen1
5050 !            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5051 !          enddo
5052 !        enddo
5053 !        write (iout,*) "fric_work and fric_work1"
5054 !        do i=1,dimen
5055 !          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5056 !        enddo
5057         do i=1,dimen
5058           do j=1,dimen
5059             ginvfric(i,j)=0.0d0
5060             do k=1,dimen
5061               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5062             enddo
5063           enddo
5064         enddo
5065         write (iout,*) "ginvfric"
5066         do i=1,dimen
5067           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5068         enddo
5069         write (iout,*) "symmetry check"
5070         do i=1,dimen
5071           do j=1,i-1
5072             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5073           enddo
5074         enddo
5075       endif
5076       return
5077       end subroutine friction_force
5078 !-----------------------------------------------------------------------------
5079       subroutine setup_fricmat
5080
5081 !     use MPI
5082       use energy_data
5083       use control_data, only:time_Bcast
5084       use control, only:tcpu
5085       use comm_syfek
5086 !      implicit real*8 (a-h,o-z)
5087 #ifdef MPI
5088       use MPI_data
5089       include 'mpif.h'
5090       real(kind=8) :: time00
5091 #endif
5092 !      include 'DIMENSIONS'
5093 !      include 'COMMON.VAR'
5094 !      include 'COMMON.CHAIN'
5095 !      include 'COMMON.DERIV'
5096 !      include 'COMMON.GEO'
5097 !      include 'COMMON.LOCAL'
5098 !      include 'COMMON.INTERACT'
5099 !      include 'COMMON.MD'
5100 !      include 'COMMON.SETUP'
5101 !      include 'COMMON.TIME1'
5102 !      integer licznik /0/
5103 !      save licznik
5104 !#ifndef LANG0
5105 !      include 'COMMON.LANGEVIN'
5106 !#else
5107 !      include 'COMMON.LANGEVIN.lang0'
5108 !#endif
5109 !      include 'COMMON.IOUNITS'
5110       integer :: IERROR
5111       integer :: i,j,ind,ind1,m
5112       logical :: lprn = .false.
5113       real(kind=8) :: dtdi !el ,gamvec(2*nres)
5114 !el      real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5115 !      real(kind=8),allocatable,dimension(:,:) :: fcopy
5116 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5117 !el      common /syfek/ gamvec
5118       real(kind=8) :: work(8*2*nres)
5119       integer :: iwork(2*nres)
5120 !el      common /przechowalnia/ ginvfric,Ghalf,fcopy
5121       integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5122       nres2=2*nres
5123       nres6=6*nres
5124 #ifdef MPI
5125       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5126        if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5127       if (fg_rank.ne.king) goto 10
5128 #endif
5129 !      nres2=2*nres
5130 !      nres6=6*nres
5131
5132       if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5133       if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5134        if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5135 !el      allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5136       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5137
5138       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5139 !  Zeroing out fricmat
5140       do i=1,dimen
5141         do j=1,dimen
5142           fricmat(i,j)=0.0d0
5143         enddo
5144       enddo
5145 !  Load the friction coefficients corresponding to peptide groups
5146       ind1=0
5147       do i=nnt,nct-1
5148         mnum=molnum(i)
5149         ind1=ind1+1
5150         gamvec(ind1)=gamp(mnum)
5151       enddo
5152 !  Load the friction coefficients corresponding to side chains
5153       m=nct-nnt
5154       ind=0
5155       do j=1,2
5156       gamsc(ntyp1_molec(j),j)=1.0d0
5157       enddo
5158       do i=nnt,nct
5159         mnum=molnum(i)
5160         ind=ind+1
5161         ii = ind+m
5162         iti=itype(i,mnum)
5163         gamvec(ii)=gamsc(iabs(iti),mnum)
5164       enddo
5165       if (surfarea) call sdarea(gamvec)
5166 !      if (lprn) then
5167 !        write (iout,*) "Matrix A and vector gamma"
5168 !        do i=1,dimen1
5169 !          write (iout,'(i2,$)') i
5170 !          do j=1,dimen
5171 !            write (iout,'(f4.1,$)') A(i,j)
5172 !          enddo
5173 !          write (iout,'(f8.3)') gamvec(i)
5174 !        enddo
5175 !      endif
5176       if (lprn) then
5177         write (iout,*) "Vector gamvec"
5178         do i=1,dimen1
5179           write (iout,'(i5,f10.5)') i, gamvec(i)
5180         enddo
5181       endif
5182
5183 ! The friction matrix
5184       do k=1,dimen
5185        do i=1,dimen
5186          dtdi=0.0d0
5187          do j=1,dimen1
5188            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5189          enddo
5190          fricmat(k,i)=dtdi
5191        enddo
5192       enddo
5193
5194       if (lprn) then
5195         write (iout,'(//a)') "Matrix fricmat"
5196         call matout2(dimen,dimen,nres2,nres2,fricmat)
5197       endif
5198       if (lang.eq.2 .or. lang.eq.3) then
5199 ! Mass-scale the friction matrix if non-direct integration will be performed
5200       do i=1,dimen
5201         do j=1,dimen
5202           Ginvfric(i,j)=0.0d0
5203           do k=1,dimen
5204             do l=1,dimen
5205               Ginvfric(i,j)=Ginvfric(i,j)+ &
5206                 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5207             enddo
5208           enddo
5209         enddo
5210       enddo
5211 ! Diagonalize the friction matrix
5212       ind=0
5213       do i=1,dimen
5214         do j=1,i
5215           ind=ind+1
5216           Ghalf(ind)=Ginvfric(i,j)
5217         enddo
5218       enddo
5219       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5220         ierr,iwork)
5221       if (lprn) then
5222         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5223           " mass-scaled friction matrix"
5224         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5225       endif
5226 ! Precompute matrices for tinker stochastic integrator
5227 #ifndef LANG0
5228       do i=1,dimen
5229         do j=1,dimen
5230           mt1(i,j)=0.0d0
5231           mt2(i,j)=0.0d0
5232           do k=1,dimen
5233             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5234             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5235           enddo
5236           mt3(j,i)=mt1(i,j)
5237         enddo
5238       enddo
5239 #endif
5240       else if (lang.eq.4) then
5241 ! Diagonalize the friction matrix
5242       ind=0
5243       do i=1,dimen
5244         do j=1,i
5245           ind=ind+1
5246           Ghalf(ind)=fricmat(i,j)
5247         enddo
5248       enddo
5249       call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5250         ierr,iwork)
5251       if (lprn) then
5252         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5253           " friction matrix"
5254         call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5255       endif
5256 ! Determine the number of zero eigenvalues of the friction matrix
5257       nzero=max0(dimen-dimen1,0)
5258 !      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5259 !        nzero=nzero+1
5260 !      enddo
5261       write (iout,*) "Number of zero eigenvalues:",nzero
5262       do i=1,dimen
5263         do j=1,dimen
5264           fricmat(i,j)=0.0d0
5265           do k=nzero+1,dimen
5266             fricmat(i,j)=fricmat(i,j) &
5267               +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5268           enddo
5269         enddo
5270       enddo
5271       if (lprn) then
5272         write (iout,'(//a)') "Generalized inverse of fricmat"
5273         call matout(dimen,dimen,nres6,nres6,fricmat)
5274       endif
5275       endif
5276 #ifdef MPI
5277   10  continue
5278       if (nfgtasks.gt.1) then
5279         if (fg_rank.eq.0) then
5280 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5281 #ifdef MPI
5282           time00=MPI_Wtime()
5283 #else
5284           time00=tcpu()
5285 #endif
5286           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5287 #ifdef MPI
5288           time_Bcast=time_Bcast+MPI_Wtime()-time00
5289 #else
5290           time_Bcast=time_Bcast+tcpu()-time00
5291 #endif
5292 !          print *,"Processor",myrank,
5293 !     &       " BROADCAST iorder in SETUP_FRICMAT"
5294         endif
5295 !       licznik=licznik+1
5296         write (iout,*) "setup_fricmat licznik"!,licznik !sp
5297 #ifdef MPI
5298         time00=MPI_Wtime()
5299 #else
5300         time00=tcpu()
5301 #endif
5302 ! Scatter the friction matrix
5303         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5304           nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5305           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5306 #ifdef TIMING
5307 #ifdef MPI
5308         time_scatter=time_scatter+MPI_Wtime()-time00
5309         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5310 #else
5311         time_scatter=time_scatter+tcpu()-time00
5312         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5313 #endif
5314 #endif
5315         do i=1,dimen
5316           do j=1,2*my_ng_count
5317             fricmat(j,i)=fcopy(i,j)
5318           enddo
5319         enddo
5320 !        write (iout,*) "My chunk of fricmat"
5321 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5322       endif
5323 #endif
5324       return
5325       end subroutine setup_fricmat
5326 !-----------------------------------------------------------------------------
5327       subroutine stochastic_force(stochforcvec)
5328
5329       use energy_data
5330       use random, only:anorm_distr
5331 !      implicit real*8 (a-h,o-z)
5332 !      include 'DIMENSIONS'
5333       use control, only: tcpu
5334       use control_data
5335 #ifdef MPI
5336       include 'mpif.h'
5337 #endif
5338 !      include 'COMMON.VAR'
5339 !      include 'COMMON.CHAIN'
5340 !      include 'COMMON.DERIV'
5341 !      include 'COMMON.GEO'
5342 !      include 'COMMON.LOCAL'
5343 !      include 'COMMON.INTERACT'
5344 !      include 'COMMON.MD'
5345 !      include 'COMMON.TIME1'
5346 !#ifndef LANG0
5347 !      include 'COMMON.LANGEVIN'
5348 !#else
5349 !      include 'COMMON.LANGEVIN.lang0'
5350 !#endif
5351 !      include 'COMMON.IOUNITS'
5352
5353       real(kind=8) :: x,sig,lowb,highb
5354       real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5355       real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5356       real(kind=8) :: time00
5357       logical :: lprn = .false.
5358       integer :: i,j,ind,mnum
5359
5360       do i=0,2*nres
5361         do j=1,3
5362           stochforc(j,i)=0.0d0
5363         enddo
5364       enddo
5365       x=0.0d0
5366
5367 #ifdef MPI
5368       time00=MPI_Wtime()
5369 #else
5370       time00=tcpu()
5371 #endif
5372 ! Compute the stochastic forces acting on bodies. Store in force.
5373       do i=nnt,nct-1
5374         sig=stdforcp(i)
5375         lowb=-5*sig
5376         highb=5*sig
5377         do j=1,3
5378           force(j,i)=anorm_distr(x,sig,lowb,highb)
5379         enddo
5380       enddo
5381       do i=nnt,nct
5382         sig2=stdforcsc(i)
5383         lowb2=-5*sig2
5384         highb2=5*sig2
5385         do j=1,3
5386           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5387         enddo
5388       enddo
5389 #ifdef MPI
5390       time_fsample=time_fsample+MPI_Wtime()-time00
5391 #else
5392       time_fsample=time_fsample+tcpu()-time00
5393 #endif
5394 ! Compute the stochastic forces acting on virtual-bond vectors.
5395       do j=1,3
5396         ff(j)=0.0d0
5397       enddo
5398       do i=nct-1,nnt,-1
5399         do j=1,3
5400           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5401         enddo
5402         do j=1,3
5403           ff(j)=ff(j)+force(j,i)
5404         enddo
5405 !        if (itype(i+1,1).ne.ntyp1) then
5406          mnum=molnum(i)
5407          if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5408           do j=1,3
5409             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5410             ff(j)=ff(j)+force(j,i+nres+1)
5411           enddo
5412         endif
5413       enddo
5414       do j=1,3
5415         stochforc(j,0)=ff(j)+force(j,nnt+nres)
5416       enddo
5417       do i=nnt,nct
5418         mnum=molnum(i)
5419         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5420         .and.(mnum.ne.5)) then
5421 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5422           do j=1,3
5423             stochforc(j,i+nres)=force(j,i+nres)
5424           enddo
5425         endif
5426       enddo
5427
5428       do j=1,3
5429         stochforcvec(j)=stochforc(j,0)
5430       enddo
5431       ind=3
5432       do i=nnt,nct-1
5433         do j=1,3
5434           stochforcvec(ind+j)=stochforc(j,i)
5435         enddo
5436         ind=ind+3
5437       enddo
5438       do i=nnt,nct
5439         mnum=molnum(i)
5440         if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5441         .and.(mnum.ne.5)) then
5442 !        if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5443           do j=1,3
5444             stochforcvec(ind+j)=stochforc(j,i+nres)
5445           enddo
5446           ind=ind+3
5447         endif
5448       enddo
5449       if (lprn) then
5450         write (iout,*) "stochforcvec"
5451         do i=1,3*dimen
5452           write(iout,'(i5,e15.5)') i,stochforcvec(i)
5453         enddo
5454         write(iout,*) "Stochastic forces backbone"
5455         do i=0,nct-1
5456           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5457         enddo
5458         write(iout,*) "Stochastic forces side chain"
5459         do i=nnt,nct
5460           write(iout,'(i5,3e15.5)') &
5461             i,(stochforc(j,i+nres),j=1,3)
5462         enddo
5463       endif
5464
5465       if (lprn) then
5466
5467       ind=0
5468       do i=nnt,nct-1
5469         write (iout,*) i,ind
5470         do j=1,3
5471           forcvec(ind+j)=force(j,i)
5472         enddo
5473         ind=ind+3
5474       enddo
5475       do i=nnt,nct
5476         write (iout,*) i,ind
5477         do j=1,3
5478           forcvec(j+ind)=force(j,i+nres)
5479         enddo
5480         ind=ind+3
5481       enddo
5482
5483       write (iout,*) "forcvec"
5484       ind=0
5485       do i=nnt,nct-1
5486         do j=1,3
5487           write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5488             forcvec(ind+j)
5489         enddo
5490         ind=ind+3
5491       enddo
5492       do i=nnt,nct
5493         do j=1,3
5494           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5495            forcvec(ind+j)
5496         enddo
5497         ind=ind+3
5498       enddo
5499
5500       endif
5501
5502       return
5503       end subroutine stochastic_force
5504 !-----------------------------------------------------------------------------
5505       subroutine sdarea(gamvec)
5506 !
5507 ! Scale the friction coefficients according to solvent accessible surface areas
5508 ! Code adapted from TINKER
5509 ! AL 9/3/04
5510 !
5511       use energy_data
5512 !      implicit real*8 (a-h,o-z)
5513 !      include 'DIMENSIONS'
5514 !      include 'COMMON.CONTROL'
5515 !      include 'COMMON.VAR'
5516 !      include 'COMMON.MD'
5517 !#ifndef LANG0
5518 !      include 'COMMON.LANGEVIN'
5519 !#else
5520 !      include 'COMMON.LANGEVIN.lang0'
5521 !#endif
5522 !      include 'COMMON.CHAIN'
5523 !      include 'COMMON.DERIV'
5524 !      include 'COMMON.GEO'
5525 !      include 'COMMON.LOCAL'
5526 !      include 'COMMON.INTERACT'
5527 !      include 'COMMON.IOUNITS'
5528 !      include 'COMMON.NAMES'
5529       real(kind=8),dimension(2*nres) :: radius,gamvec   !(maxres2)
5530       real(kind=8),parameter :: twosix = 1.122462048309372981d0
5531       logical :: lprn = .false.
5532       real(kind=8) :: probe,area,ratio
5533       integer :: i,j,ind,iti,mnum
5534 !
5535 !     determine new friction coefficients every few SD steps
5536 !
5537 !     set the atomic radii to estimates of sigma values
5538 !
5539 !      print *,"Entered sdarea"
5540       probe = 0.0d0
5541
5542       do i=1,2*nres
5543         radius(i)=0.0d0
5544       enddo
5545 !  Load peptide group radii
5546       do i=nnt,nct-1
5547         mnum=molnum(i)
5548         radius(i)=pstok(mnum)
5549       enddo
5550 !  Load side chain radii
5551       do i=nnt,nct
5552         mnum=molnum(i)
5553         iti=itype(i,mnum)
5554         radius(i+nres)=restok(iti,mnum)
5555       enddo
5556 !      do i=1,2*nres
5557 !        write (iout,*) "i",i," radius",radius(i)
5558 !      enddo
5559       do i = 1, 2*nres
5560          radius(i) = radius(i) / twosix
5561          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
5562       end do
5563 !
5564 !     scale atomic friction coefficients by accessible area
5565 !
5566       if (lprn) write (iout,*) &
5567         "Original gammas, surface areas, scaling factors, new gammas, ",&
5568         "std's of stochastic forces"
5569       ind=0
5570       do i=nnt,nct-1
5571         if (radius(i).gt.0.0d0) then
5572           call surfatom (i,area,radius)
5573           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5574           if (lprn) write (iout,'(i5,3f10.5,$)') &
5575             i,gamvec(ind+1),area,ratio
5576           do j=1,3
5577             ind=ind+1
5578             gamvec(ind) = ratio * gamvec(ind)
5579           enddo
5580           mnum=molnum(i)
5581           stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5582           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5583         endif
5584       enddo
5585       do i=nnt,nct
5586         if (radius(i+nres).gt.0.0d0) then
5587           call surfatom (i+nres,area,radius)
5588           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5589           if (lprn) write (iout,'(i5,3f10.5,$)') &
5590             i,gamvec(ind+1),area,ratio
5591           do j=1,3
5592             ind=ind+1
5593             gamvec(ind) = ratio * gamvec(ind)
5594           enddo
5595           mnum=molnum(i)
5596           stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5597           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5598         endif
5599       enddo
5600
5601       return
5602       end subroutine sdarea
5603 !-----------------------------------------------------------------------------
5604 ! surfatom.f
5605 !-----------------------------------------------------------------------------
5606 !
5607 !
5608 !     ###################################################
5609 !     ##  COPYRIGHT (C)  1996  by  Jay William Ponder  ##
5610 !     ##              All Rights Reserved              ##
5611 !     ###################################################
5612 !
5613 !     ################################################################
5614 !     ##                                                            ##
5615 !     ##  subroutine surfatom  --  exposed surface area of an atom  ##
5616 !     ##                                                            ##
5617 !     ################################################################
5618 !
5619 !
5620 !     "surfatom" performs an analytical computation of the surface
5621 !     area of a specified atom; a simplified version of "surface"
5622 !
5623 !     literature references:
5624 !
5625 !     T. J. Richmond, "Solvent Accessible Surface Area and
5626 !     Excluded Volume in Proteins", Journal of Molecular Biology,
5627 !     178, 63-89 (1984)
5628 !
5629 !     L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5630 !     Applied to Molecular Dynamics of Proteins in Solution",
5631 !     Protein Science, 1, 227-235 (1992)
5632 !
5633 !     variables and parameters:
5634 !
5635 !     ir       number of atom for which area is desired
5636 !     area     accessible surface area of the atom
5637 !     radius   radii of each of the individual atoms
5638 !
5639 !
5640       subroutine surfatom(ir,area,radius)
5641
5642 !      implicit real*8 (a-h,o-z)
5643 !      include 'DIMENSIONS'
5644 !      include 'sizes.i'
5645 !      include 'COMMON.GEO'
5646 !      include 'COMMON.IOUNITS'
5647 !      integer :: nres,
5648       integer :: nsup,nstart_sup
5649 !      double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5650 !      common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5651 !     & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5652 !     & dc_work(MAXRES6),nres,nres0
5653       integer,parameter :: maxarc=300
5654       integer :: i,j,k,m
5655       integer :: ii,ib,jb
5656       integer :: io,ir
5657       integer :: mi,ni,narc
5658       integer :: key(maxarc)
5659       integer :: intag(maxarc)
5660       integer :: intag1(maxarc)
5661       real(kind=8) :: area,arcsum
5662       real(kind=8) :: arclen,exang
5663       real(kind=8) :: delta,delta2
5664       real(kind=8) :: eps,rmove
5665       real(kind=8) :: xr,yr,zr
5666       real(kind=8) :: rr,rrsq
5667       real(kind=8) :: rplus,rminus
5668       real(kind=8) :: axx,axy,axz
5669       real(kind=8) :: ayx,ayy
5670       real(kind=8) :: azx,azy,azz
5671       real(kind=8) :: uxj,uyj,uzj
5672       real(kind=8) :: tx,ty,tz
5673       real(kind=8) :: txb,tyb,td
5674       real(kind=8) :: tr2,tr,txr,tyr
5675       real(kind=8) :: tk1,tk2
5676       real(kind=8) :: thec,the,t,tb
5677       real(kind=8) :: txk,tyk,tzk
5678       real(kind=8) :: t1,ti,tf,tt
5679       real(kind=8) :: txj,tyj,tzj
5680       real(kind=8) :: ccsq,cc,xysq
5681       real(kind=8) :: bsqk,bk,cosine
5682       real(kind=8) :: dsqj,gi,pix2
5683       real(kind=8) :: therk,dk,gk
5684       real(kind=8) :: risqk,rik
5685       real(kind=8) :: radius(2*nres)    !(maxatm) (maxatm=maxres2)
5686       real(kind=8) :: ri(maxarc),risq(maxarc)
5687       real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5688       real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5689       real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5690       real(kind=8) :: dsq(maxarc),bsq(maxarc)
5691       real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5692       real(kind=8) :: arci(maxarc),arcf(maxarc)
5693       real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5694       real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5695       real(kind=8) :: kent(maxarc),kout(maxarc)
5696       real(kind=8) :: ther(maxarc)
5697       logical :: moved,top
5698       logical :: omit(maxarc)
5699 !
5700 !      include 'sizes.i'
5701 !      maxatm = 2*nres  !maxres2 maxres2=2*maxres
5702 !      maxlight = 8*maxatm
5703 !      maxbnd = 2*maxatm
5704 !      maxang = 3*maxatm
5705 !      maxtors = 4*maxatm
5706 !
5707 !     zero out the surface area for the sphere of interest
5708 !
5709       area = 0.0d0
5710 !      write (2,*) "ir",ir," radius",radius(ir)
5711       if (radius(ir) .eq. 0.0d0)  return
5712 !
5713 !     set the overlap significance and connectivity shift
5714 !
5715       pix2 = 2.0d0 * pi
5716       delta = 1.0d-8
5717       delta2 = delta * delta
5718       eps = 1.0d-8
5719       moved = .false.
5720       rmove = 1.0d-8
5721 !
5722 !     store coordinates and radius of the sphere of interest
5723 !
5724       xr = c(1,ir)
5725       yr = c(2,ir)
5726       zr = c(3,ir)
5727       rr = radius(ir)
5728       rrsq = rr * rr
5729 !
5730 !     initialize values of some counters and summations
5731 !
5732    10 continue
5733       io = 0
5734       jb = 0
5735       ib = 0
5736       arclen = 0.0d0
5737       exang = 0.0d0
5738 !
5739 !     test each sphere to see if it overlaps the sphere of interest
5740 !
5741       do i = 1, 2*nres
5742          if (i.eq.ir .or. radius(i).eq.0.0d0)  goto 30
5743          rplus = rr + radius(i)
5744          tx = c(1,i) - xr
5745          if (abs(tx) .ge. rplus)  goto 30
5746          ty = c(2,i) - yr
5747          if (abs(ty) .ge. rplus)  goto 30
5748          tz = c(3,i) - zr
5749          if (abs(tz) .ge. rplus)  goto 30
5750 !
5751 !     check for sphere overlap by testing distance against radii
5752 !
5753          xysq = tx*tx + ty*ty
5754          if (xysq .lt. delta2) then
5755             tx = delta
5756             ty = 0.0d0
5757             xysq = delta2
5758          end if
5759          ccsq = xysq + tz*tz
5760          cc = sqrt(ccsq)
5761          if (rplus-cc .le. delta)  goto 30
5762          rminus = rr - radius(i)
5763 !
5764 !     check to see if sphere of interest is completely buried
5765 !
5766          if (cc-abs(rminus) .le. delta) then
5767             if (rminus .le. 0.0d0)  goto 170
5768             goto 30
5769          end if
5770 !
5771 !     check for too many overlaps with sphere of interest
5772 !
5773          if (io .ge. maxarc) then
5774             write (iout,20)
5775    20       format (/,' SURFATOM  --  Increase the Value of MAXARC')
5776             stop
5777          end if
5778 !
5779 !     get overlap between current sphere and sphere of interest
5780 !
5781          io = io + 1
5782          xc1(io) = tx
5783          yc1(io) = ty
5784          zc1(io) = tz
5785          dsq1(io) = xysq
5786          bsq1(io) = ccsq
5787          b1(io) = cc
5788          gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5789          intag1(io) = i
5790          omit(io) = .false.
5791    30    continue
5792       end do
5793 !
5794 !     case where no other spheres overlap the sphere of interest
5795 !
5796       if (io .eq. 0) then
5797          area = 4.0d0 * pi * rrsq
5798          return
5799       end if
5800 !
5801 !     case where only one sphere overlaps the sphere of interest
5802 !
5803       if (io .eq. 1) then
5804          area = pix2 * (1.0d0 + gr(1))
5805          area = mod(area,4.0d0*pi) * rrsq
5806          return
5807       end if
5808 !
5809 !     case where many spheres intersect the sphere of interest;
5810 !     sort the intersecting spheres by their degree of overlap
5811 !
5812       call sort2 (io,gr,key)
5813       do i = 1, io
5814          k = key(i)
5815          intag(i) = intag1(k)
5816          xc(i) = xc1(k)
5817          yc(i) = yc1(k)
5818          zc(i) = zc1(k)
5819          dsq(i) = dsq1(k)
5820          b(i) = b1(k)
5821          bsq(i) = bsq1(k)
5822       end do
5823 !
5824 !     get radius of each overlap circle on surface of the sphere
5825 !
5826       do i = 1, io
5827          gi = gr(i) * rr
5828          bg(i) = b(i) * gi
5829          risq(i) = rrsq - gi*gi
5830          ri(i) = sqrt(risq(i))
5831          ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5832       end do
5833 !
5834 !     find boundary of inaccessible area on sphere of interest
5835 !
5836       do k = 1, io-1
5837          if (.not. omit(k)) then
5838             txk = xc(k)
5839             tyk = yc(k)
5840             tzk = zc(k)
5841             bk = b(k)
5842             therk = ther(k)
5843 !
5844 !     check to see if J circle is intersecting K circle;
5845 !     get distance between circle centers and sum of radii
5846 !
5847             do j = k+1, io
5848                if (omit(j))  goto 60
5849                cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5850                cc = acos(min(1.0d0,max(-1.0d0,cc)))
5851                td = therk + ther(j)
5852 !
5853 !     check to see if circles enclose separate regions
5854 !
5855                if (cc .ge. td)  goto 60
5856 !
5857 !     check for circle J completely inside circle K
5858 !
5859                if (cc+ther(j) .lt. therk)  goto 40
5860 !
5861 !     check for circles that are essentially parallel
5862 !
5863                if (cc .gt. delta)  goto 50
5864    40          continue
5865                omit(j) = .true.
5866                goto 60
5867 !
5868 !     check to see if sphere of interest is completely buried
5869 !
5870    50          continue
5871                if (pix2-cc .le. td)  goto 170
5872    60          continue
5873             end do
5874          end if
5875       end do
5876 !
5877 !     find T value of circle intersections
5878 !
5879       do k = 1, io
5880          if (omit(k))  goto 110
5881          omit(k) = .true.
5882          narc = 0
5883          top = .false.
5884          txk = xc(k)
5885          tyk = yc(k)
5886          tzk = zc(k)
5887          dk = sqrt(dsq(k))
5888          bsqk = bsq(k)
5889          bk = b(k)
5890          gk = gr(k) * rr
5891          risqk = risq(k)
5892          rik = ri(k)
5893          therk = ther(k)
5894 !
5895 !     rotation matrix elements
5896 !
5897          t1 = tzk / (bk*dk)
5898          axx = txk * t1
5899          axy = tyk * t1
5900          axz = dk / bk
5901          ayx = tyk / dk
5902          ayy = txk / dk
5903          azx = txk / bk
5904          azy = tyk / bk
5905          azz = tzk / bk
5906          do j = 1, io
5907             if (.not. omit(j)) then
5908                txj = xc(j)
5909                tyj = yc(j)
5910                tzj = zc(j)
5911 !
5912 !     rotate spheres so K vector colinear with z-axis
5913 !
5914                uxj = txj*axx + tyj*axy - tzj*axz
5915                uyj = tyj*ayy - txj*ayx
5916                uzj = txj*azx + tyj*azy + tzj*azz
5917                cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5918                if (acos(cosine) .lt. therk+ther(j)) then
5919                   dsqj = uxj*uxj + uyj*uyj
5920                   tb = uzj*gk - bg(j)
5921                   txb = uxj * tb
5922                   tyb = uyj * tb
5923                   td = rik * dsqj
5924                   tr2 = risqk*dsqj - tb*tb
5925                   tr2 = max(eps,tr2)
5926                   tr = sqrt(tr2)
5927                   txr = uxj * tr
5928                   tyr = uyj * tr
5929 !
5930 !     get T values of intersection for K circle
5931 !
5932                   tb = (txb+tyr) / td
5933                   tb = min(1.0d0,max(-1.0d0,tb))
5934                   tk1 = acos(tb)
5935                   if (tyb-txr .lt. 0.0d0)  tk1 = pix2 - tk1
5936                   tb = (txb-tyr) / td
5937                   tb = min(1.0d0,max(-1.0d0,tb))
5938                   tk2 = acos(tb)
5939                   if (tyb+txr .lt. 0.0d0)  tk2 = pix2 - tk2
5940                   thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5941                   if (abs(thec) .lt. 1.0d0) then
5942                      the = -acos(thec)
5943                   else if (thec .ge. 1.0d0) then
5944                      the = 0.0d0
5945                   else if (thec .le. -1.0d0) then
5946                      the = -pi
5947                   end if
5948 !
5949 !     see if "tk1" is entry or exit point; check t=0 point;
5950 !     "ti" is exit point, "tf" is entry point
5951 !
5952                   cosine = min(1.0d0,max(-1.0d0, &
5953                                   (uzj*gk-uxj*rik)/(b(j)*rr)))
5954                   if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5955                      ti = tk2
5956                      tf = tk1
5957                   else
5958                      ti = tk2
5959                      tf = tk1
5960                   end if
5961                   narc = narc + 1
5962                   if (narc .ge. maxarc) then
5963                      write (iout,70)
5964    70                format (/,' SURFATOM  --  Increase the Value',&
5965                                 ' of MAXARC')
5966                      stop
5967                   end if
5968                   if (tf .le. ti) then
5969                      arcf(narc) = tf
5970                      arci(narc) = 0.0d0
5971                      tf = pix2
5972                      lt(narc) = j
5973                      ex(narc) = the
5974                      top = .true.
5975                      narc = narc + 1
5976                   end if
5977                   arcf(narc) = tf
5978                   arci(narc) = ti
5979                   lt(narc) = j
5980                   ex(narc) = the
5981                   ux(j) = uxj
5982                   uy(j) = uyj
5983                   uz(j) = uzj
5984                end if
5985             end if
5986          end do
5987          omit(k) = .false.
5988 !
5989 !     special case; K circle without intersections
5990 !
5991          if (narc .le. 0)  goto 90
5992 !
5993 !     general case; sum up arclength and set connectivity code
5994 !
5995          call sort2 (narc,arci,key)
5996          arcsum = arci(1)
5997          mi = key(1)
5998          t = arcf(mi)
5999          ni = mi
6000          if (narc .gt. 1) then
6001             do j = 2, narc
6002                m = key(j)
6003                if (t .lt. arci(j)) then
6004                   arcsum = arcsum + arci(j) - t
6005                   exang = exang + ex(ni)
6006                   jb = jb + 1
6007                   if (jb .ge. maxarc) then
6008                      write (iout,80)
6009    80                format (/,' SURFATOM  --  Increase the Value',&
6010                                 ' of MAXARC')
6011                      stop
6012                   end if
6013                   i = lt(ni)
6014                   kent(jb) = maxarc*i + k
6015                   i = lt(m)
6016                   kout(jb) = maxarc*k + i
6017                end if
6018                tt = arcf(m)
6019                if (tt .ge. t) then
6020                   t = tt
6021                   ni = m
6022                end if
6023             end do
6024          end if
6025          arcsum = arcsum + pix2 - t
6026          if (.not. top) then
6027             exang = exang + ex(ni)
6028             jb = jb + 1
6029             i = lt(ni)
6030             kent(jb) = maxarc*i + k
6031             i = lt(mi)
6032             kout(jb) = maxarc*k + i
6033          end if
6034          goto 100
6035    90    continue
6036          arcsum = pix2
6037          ib = ib + 1
6038   100    continue
6039          arclen = arclen + gr(k)*arcsum
6040   110    continue
6041       end do
6042       if (arclen .eq. 0.0d0)  goto 170
6043       if (jb .eq. 0)  goto 150
6044 !
6045 !     find number of independent boundaries and check connectivity
6046 !
6047       j = 0
6048       do k = 1, jb
6049          if (kout(k) .ne. 0) then
6050             i = k
6051   120       continue
6052             m = kout(i)
6053             kout(i) = 0
6054             j = j + 1
6055             do ii = 1, jb
6056                if (m .eq. kent(ii)) then
6057                   if (ii .eq. k) then
6058                      ib = ib + 1
6059                      if (j .eq. jb)  goto 150
6060                      goto 130
6061                   end if
6062                   i = ii
6063                   goto 120
6064                end if
6065             end do
6066   130       continue
6067          end if
6068       end do
6069       ib = ib + 1
6070 !
6071 !     attempt to fix connectivity error by moving atom slightly
6072 !
6073       if (moved) then
6074          write (iout,140)  ir
6075   140    format (/,' SURFATOM  --  Connectivity Error at Atom',i6)
6076       else
6077          moved = .true.
6078          xr = xr + rmove
6079          yr = yr + rmove
6080          zr = zr + rmove
6081          goto 10
6082       end if
6083 !
6084 !     compute the exposed surface area for the sphere of interest
6085 !
6086   150 continue
6087       area = ib*pix2 + exang + arclen
6088       area = mod(area,4.0d0*pi) * rrsq
6089 !
6090 !     attempt to fix negative area by moving atom slightly
6091 !
6092       if (area .lt. 0.0d0) then
6093          if (moved) then
6094             write (iout,160)  ir
6095   160       format (/,' SURFATOM  --  Negative Area at Atom',i6)
6096          else
6097             moved = .true.
6098             xr = xr + rmove
6099             yr = yr + rmove
6100             zr = zr + rmove
6101             goto 10
6102          end if
6103       end if
6104   170 continue
6105       return
6106       end subroutine surfatom
6107 !----------------------------------------------------------------
6108 !----------------------------------------------------------------
6109       subroutine alloc_MD_arrays
6110 !EL Allocation of arrays used by MD module
6111
6112       integer :: nres2,nres6
6113       nres2=nres*2
6114       nres6=nres*6
6115 !----------------------
6116 #ifndef LANG0
6117 ! commom.langevin
6118 !      common /langforc/
6119       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6120       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6121       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6122       allocate(fricvec(nres2,nres2))
6123       allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6124       allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6125       allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6126       allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6127       allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6128       allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6129       allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6130       allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6131       allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6132       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6133 !      common /langmat/
6134       allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6135 !----------------------
6136 #else
6137 ! commom.langevin.lang0
6138 !      common /langforc/
6139       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6140       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6141       allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6142       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6143       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6144 #endif
6145
6146       if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6147 !----------------------
6148 ! commom.hairpin in CSA module
6149 !----------------------
6150 ! common.mce in MCM_MD module
6151 !----------------------
6152 ! common.MD
6153 !      common /mdgrad/ in module.energy
6154 !      common /back_constr/ in module.energy
6155 !      common /qmeas/ in module.energy
6156 !      common /mdpar/
6157 !      common /MDcalc/
6158       allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6159 !      common /lagrange/
6160       allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6161       allocate(d_a_work(nres6)) !(6*MAXRES)
6162 #ifdef FIVEDIAG
6163       allocate(DM(nres2),DU1(nres2),DU2(nres2))
6164       allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6165 #else
6166       write (iout,*) "Before A Allocation",nfgtasks-1
6167       call flush(iout)
6168       allocate(Gmat(nres2,nres2),A(nres2,nres2))
6169       if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6170       allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6171 #endif
6172       allocate(Geigen(nres2)) !(maxres2)
6173       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6174 !      common /inertia/ in io_conf: parmread
6175 !      real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6176 !      common /langevin/in io read_MDpar
6177 !      real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6178 !      real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6179 ! in io_conf: parmread
6180 !      real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6181 !      common /mdpmpi/ in control: ergastulum
6182       if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6183       if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6184       if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6185       if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6186 !----------------------
6187 ! common.muca in read_muca
6188 !      common /double_muca/
6189 !      real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6190 !      real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6191 !       nemuca2,hist !(4*maxres)
6192 !      common /integer_muca/
6193 !      integer :: nmuca,imtime,muca_smooth
6194 !      common /mucarem/
6195 !      real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6196 !----------------------
6197 ! common.MD
6198 !      common /mdgrad/ in module.energy
6199 !      common /back_constr/ in module.energy
6200 !      common /qmeas/ in module.energy
6201 !      common /mdpar/
6202 !      common /MDcalc/
6203 !      common /lagrange/
6204       allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6205       allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6206       allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6207       allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6208 !----------------------
6209 !      COMMON /BANII/ D
6210       allocate(D_ban(nres6))    !(MAXRES6) maxres6=6*maxres
6211 !      common /stochcalc/ stochforcvec
6212       allocate(stochforcvec(nres6))     !(MAXRES6) maxres6=6*maxres
6213 !----------------------
6214       return
6215       end subroutine alloc_MD_arrays
6216 !-----------------------------------------------------------------------------
6217 !-----------------------------------------------------------------------------
6218       end module MDyn