source/unres/MD.f90

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