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