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