source/unres/MD.F90

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