source/unres/MD.f90

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