source/unres/MD.f90

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