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