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