source/unres/MD.f90

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