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