source/unres/MD.f90

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