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