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