source/unres/MD.F90

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