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