source/unres/MD.f90

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