source/unres/REMD.F90

   1       module REMD
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use MD_data
   5       use REMD_data
   6       use muca_md
   7
   8       implicit none
   9 !-----------------------------------------------------------------------------
  10 !
  11 !
  12 !-----------------------------------------------------------------------------
  13       contains
  14 !-----------------------------------------------------------------------------
  15 ! lagrangian_lesyng.F
  16 !-----------------------------------------------------------------------------
  17       subroutine lagrangian
  18 !-------------------------------------------------------------------------
  19 !  This subroutine contains the total lagrangain from which the accelerations
  20 !  are obtained.  For numerical gradient checking, the derivetive of the
  21 !  lagrangian in the velocities and coordinates are calculated seperately
  22 !-------------------------------------------------------------------------
  23 !       implicit real*8 (a-h,o-z)
  24 !       include 'DIMENSIONS'
  25       use comm_cipiszcze
  26       use energy_data
  27 #ifdef FIVEDIAG
  28       use energy, only: grad_transform
  29 #endif
  30       use geometry_data, only: nres
  31       use control_data    !el, only: mucadyn,lmuca
  32 #ifdef MPI
  33        include 'mpif.h'
  34       real(kind=8) :: time00
  35 #endif
  36 !       include 'COMMON.VAR'
  37 !       include 'COMMON.CHAIN'
  38 !       include 'COMMON.DERIV'
  39 !       include 'COMMON.GEO'
  40 !       include 'COMMON.LOCAL'
  41 !       include 'COMMON.INTERACT'
  42 !       include 'COMMON.MD'
  43 !       include 'COMMON.IOUNITS'
  44 !       include 'COMMON.CONTROL'
  45 !       include 'COMMON.MUCA'
  46 !      include 'COMMON.TIME1'
  47        integer ::i,j,ind,itime,mnum,innt,inct,inct_prev,ichain,n,mark
  48        real(kind=8) :: zapas(6*nres) !,muca_factor      !maxres6=6*maxres
  49 ! the line below might be wrong
  50 #ifdef FIVEDIAG
  51        real(kind=8) :: rs(2*nres),xsolv(2*nres)
  52 !#ifdef CHECK5DSOL
  53        real(kind=8) :: rscheck(2*nres),rsold(2*nres)
  54 !#endif
  55 #endif
  56        logical :: lprn = .false.
  57 !el       common /cipiszcze/ itime
  58        itime = itt_comm
  59
  60 #ifdef TIMING
  61        time00=MPI_Wtime()
  62 #endif
  63 #ifdef FIVEDIAG
  64       call grad_transform
  65       d_a=0.0d0
  66       if (lprn) then
  67         write (iout,*) "Potential forces backbone"
  68         do i=1,nres
  69           write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
  70         enddo
  71         write (iout,*) "Potential forces sidechain"
  72         do i=nnt,nct
  73 !          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
  74           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
  75         enddo
  76       endif
  77       do ichain=1,nchain
  78         n=dimen_chain(ichain)
  79         innt=iposd_chain(ichain)
  80         do j=1,3
  81           ind=1
  82           do i=chain_border(1,ichain),chain_border(2,ichain)
  83             mnum=molnum(i)
  84             if (itype(i,1).eq.10.or.mnum.ge.3)then
  85               rs(ind)=-gcart(j,i)-gxcart(j,i)
  86               ind=ind+1
  87             else
  88               rs(ind)=-gcart(j,i)
  89               rs(ind+1)=-gxcart(j,i)
  90               ind=ind+2
  91             end if
  92           enddo
  93
  94 #ifdef CHECK5DSOL
  95          rsold=rs
  96 #endif
  97          if (lprn) then
  98            write(iout,*) "RHS of the 5-diag equations system",&
  99            ichain," j",j
 100            do i=1,n
 101              write(iout,*) i,rs(i)
 102            enddo
 103          endif
 104
 105          call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv)
 106          if (lprn) then
 107            write (iout,*) "Solution of the 5-diagonal equations system"
 108            do i=1,n
 109              write (iout,'(i5,f10.5)') i,xsolv(i)
 110            enddo
 111          endif
 112 #ifdef CHECK5DSOL
 113 ! Check the solution
 114           call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),&
 115            xsolv,rscheck)
 116           do i=1,n
 117             write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
 118             "ratio",rscheck(i)/rsold(i)
 119           enddo
 120 ! end check
 121 #endif
 122 #undef CHECK5DSOL
 123           ind=1
 124           do i=chain_border(1,ichain),chain_border(2,ichain)
 125             mnum=molnum(i)
 126             if (itype(i,1).eq.10 .or.mnum.ge.3) then
 127               d_a(j,i)=xsolv(ind)
 128               ind=ind+1
 129             else
 130               d_a(j,i)=xsolv(ind)
 131               d_a(j,i+nres)=xsolv(ind+1)
 132               ind=ind+2
 133             end if
 134           enddo
 135         enddo ! j
 136       enddo ! ichain
 137       if (lprn) then
 138         write (iout,*) "Acceleration in CA and SC oordinates"
 139         do i=1,nres
 140           write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
 141         enddo
 142         do i=nnt,nct
 143           write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
 144         enddo
 145       endif
 146 !C Conevert d_a to virtual-bon-vector basis
 147 #define WLOS
 148 #ifdef WLOS
 149 !c      write (iout,*) "WLOS"
 150       if (nnt.eq.1) then
 151         d_a(:,0)=d_a(:,1)
 152       endif
 153       do i=1,nres
 154         mnum=molnum(i)
 155         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum) .or.mnum.ge.3) then
 156           do j=1,3
 157             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 158           enddo
 159         else
 160           do j=1,3
 161             d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 162             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 163           enddo
 164         end if
 165       enddo
 166       d_a(:,nres)=0.0d0
 167       d_a(:,nct)=0.0d0
 168       d_a(:,2*nres)=0.0d0
 169 !c      d_a(:,0)=d_a(:,1)
 170 !c      d_a(:,1)=0.0d0
 171 !c      write (iout,*) "Shifting accelerations"
 172       if (nnt.gt.1) then
 173         d_a(:,0)=d_a(:,1)
 174         d_a(:,1)=0.0d0
 175       endif
 176 #define CHUJ
 177 #ifdef CHUJ
 178       do ichain=2,nchain
 179 !TEST 27.06.2023 godz 16.00
 180       if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
 181 !c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 182 !c     &     chain_border1(1,ichain)
 183         d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
 184         d_a(:,chain_border1(1,ichain))=0.0d0
 185       enddo
 186 !c      write (iout,*) "Adding accelerations"
 187       do ichain=2,nchain
 188       if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
 189 !c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
 190 !c     &   chain_border(2,ichain-1)
 191         d_a(:,chain_border1(1,ichain)-1)=&
 192        d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
 193         d_a(:,chain_border(2,ichain-1))=0.0d0
 194       enddo
 195 #endif
 196 #else
 197       inct_prev=0
 198       do j=1,3
 199         aaux(j)=0.0d0
 200       enddo
 201       do ichain=1,nchain
 202         innt=chain_border(1,ichain)
 203         inct=chain_border(2,ichain)
 204         do j=1,3
 205           d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
 206         enddo
 207         inct_prev=inct+1
 208         do i=innt,inct
 209           mnum=molnum(i)
 210           if ((itype(i).ne.10).and.(mnum.le.3)) then
 211             do j=1,3
 212               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 213             enddo
 214           endif
 215         enddo
 216         do j=1,3
 217           aaux(j)=d_a(j,inct)
 218         enddo
 219         do i=innt,inct
 220           do j=1,3
 221             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 222           enddo
 223         enddo
 224       enddo
 225 #endif
 226       if (lprn) then
 227         write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
 228         do i=0,nres
 229           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 230         enddo
 231         do i=nnt,nct
 232           write (iout,'(i3,3f10.5,3x,3f10.5)')&
 233           i,(d_a(j,i+nres),j=1,3)
 234         enddo
 235       endif
 236
 237 #else
 238 ! Old procedure
 239        do j=1,3
 240          zapas(j)=-gcart(j,0)
 241        enddo
 242       ind=3
 243       if (lprn) then
 244         write (iout,*) "Potential forces backbone"
 245       endif
 246       do i=nnt,nct-1
 247         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 248           i,(-gcart(j,i),j=1,3)
 249         do j=1,3
 250           ind=ind+1
 251           zapas(ind)=-gcart(j,i)
 252         enddo
 253       enddo
 254       if (lprn) write (iout,*) "Potential forces sidechain"
 255       do i=nnt,nct
 256           mnum=molnum(i)
 257         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
 258           .and.mnum.lt.4) then
 259           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
 260              i,(-gxcart(j,i),j=1,3)
 261           do j=1,3
 262             ind=ind+1
 263             zapas(ind)=-gxcart(j,i)
 264           enddo
 265         endif
 266       enddo
 267
 268       call ginv_mult(zapas,d_a_work)
 269
 270       do j=1,3
 271         d_a(j,0)=d_a_work(j)
 272       enddo
 273       ind=3
 274       do i=nnt,nct-1
 275         do j=1,3
 276           ind=ind+1
 277           d_a(j,i)=d_a_work(ind)
 278         enddo
 279       enddo
 280       do i=nnt,nct
 281        mnum=molnum(i)
 282 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
 283          if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
 284           .and.mnum.lt.4) then
 285           do j=1,3
 286             ind=ind+1
 287             d_a(j,i+nres)=d_a_work(ind)
 288           enddo
 289         endif
 290       enddo
 291
 292 #endif
 293       if(lmuca) then
 294        imtime=imtime+1
 295        if(mucadyn.gt.0) call muca_update(potE)
 296        factor=muca_factor(potE)*t_bath*Rb
 297
 298 !d       print *,'lmuca ',factor,potE
 299        do j=1,3
 300           d_a(j,0)=d_a(j,0)*factor
 301        enddo
 302        do i=nnt,nct-1
 303          do j=1,3
 304           d_a(j,i)=d_a(j,i)*factor
 305          enddo
 306        enddo
 307        do i=nnt,nct
 308          do j=1,3
 309           d_a(j,i+nres)=d_a(j,i+nres)*factor
 310          enddo
 311        enddo
 312
 313       endif
 314
 315       if (lprn) then
 316         write(iout,*) 'acceleration 3D'
 317         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 318         do i=nnt,nct-1
 319          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 320         enddo
 321         do i=nnt,nct
 322          write (iout,'(i3,3f10.5,3x,3f10.5)') &
 323            i+nres,(d_a(j,i+nres),j=1,3)
 324         enddo
 325       endif
 326 #ifdef TIMING
 327       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 328 #endif
 329       return
 330       end subroutine lagrangian
 331 !-----------------------------------------------------------------------------
 332       subroutine setup_MD_matrices
 333
 334       use geometry_data, only: nres,nside
 335       use control_data
 336       use MPI_data
 337       use energy_data
 338       use geometry, only:int_bounds
 339       use md_calc
 340 !      implicit real*8 (a-h,o-z)
 341 !      include 'DIMENSIONS'
 342 #ifdef MPI
 343       include 'mpif.h'
 344       integer :: ierror
 345       real(kind=8) :: time00
 346 !      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 347 #endif
 348 !      include 'COMMON.SETUP'
 349 !      include 'COMMON.VAR'
 350 !      include 'COMMON.CHAIN'
 351 !      include 'COMMON.DERIV'
 352 !      include 'COMMON.GEO'
 353 !      include 'COMMON.LOCAL'
 354 !      include 'COMMON.INTERACT'
 355 !      include 'COMMON.MD'
 356 !#ifndef LANG0
 357 !      include 'COMMON.LANGEVIN'
 358 !#else
 359 !      include 'COMMON.LANGEVIN.lang0'
 360 !#endif
 361 !      include 'COMMON.IOUNITS'
 362 !      include 'COMMON.TIME1'
 363       logical :: lprn = .false.
 364       logical :: osob
 365 #ifndef FIVEDIAG
 366       real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
 367 #endif
 368       real(kind=8) :: dtdi
 369       real(kind=8),dimension(:),allocatable :: massvec,sqreig   !(maxres2) maxres2=2*maxres
 370       real(kind=8) :: relfeh,eps1,eps2
 371 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 372 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 373 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
 374 !el      real(kind=8),dimension(:,:),allocatable :: Gcopy
 375       real(kind=8),dimension(:),allocatable :: work     !(8*maxres6)
 376       integer,dimension(:),allocatable :: iwork !(maxres6) maxres6=6*maxres
 377 !      common /jakistam/ iwork,work,massvec,sqreig
 378 !el      common /przechowalnia/ Gcopy,Ghalf
 379       real(kind=8) :: coeff,mscab
 380       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
 381       real(kind=8) :: ip4
 382       integer :: iz,mnum,ichain,n,dimenp,innt,inct
 383       if(.not.allocated(massvec)) then
 384       allocate(massvec(2*nres),sqreig(2*nres))
 385       allocate(work(8*6*nres))
 386       allocate(iwork(6*nres))
 387       endif
 388       print *,"just entered"
 389       relfeh=1.0d-14
 390       nres2=2*nres
 391       print *,"FIVEDIAG"
 392       write (iout,*) "before FIVEDIAG"
 393 #ifdef FIVEDIAG
 394       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 395 #endif
 396 #ifndef FIVEDIAG
 397       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 398       write (iout,*) "ALLOCATE"
 399       print *,"ALLOCATE"
 400       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
 401 !      if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 402       if(.not.allocated(Bmat))  allocate(Bmat(nres2,nres2))
 403       if(.not.allocated(matmult)) allocate(matmult(nres2,nres2))
 404 #endif
 405 !
 406 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 407 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 408 !
 409 ! Determine the number of degrees of freedom (dimen) and the number of
 410 ! sites (dimen1)
 411 !      dimen=(nct-nnt+1)+nside
 412 !      dimen1=(nct-nnt)+(nct-nnt+1)
 413 !      dimen3=dimen*3
 414 !      write (iout,*) "nnt",nnt," nct",nct," nside",nside
 415 #ifdef FIVEDIAG
 416 #ifdef CRYST_BOND
 417        ip4=ip/4
 418        DM(1)=mp/4+ip4
 419        if (iabs(itype(nnt)).eq.10) then
 420          DM(1)=DM(1)+msc(10)
 421          ind=2
 422          nind=1
 423        else
 424          DM(1)=DM(1)+isc(iabs(itype(nnt)))
 425          DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
 426          ind=3
 427          nind=2
 428        endif
 429        do i=nnt+1,nct-1
 430         mnum=molnum(i)
 431 !         if (iabs(itype(i,1)).eq.ntyp1) cycle
 432          DM(ind)=2*ip4+mp/2
 433         if (iabs(itype(i,1)).eq.10 .or. &
 434           iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then
 435            if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
 436            ind=ind+1
 437          else
 438            DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
 439            DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
 440            ind=ind+2
 441          endif
 442        enddo
 443        if (nct.gt.nnt) then
 444        DM(ind)=ip4+mp/4
 445        if (iabs(itype(nct)).eq.10) then
 446          DM(ind)=DM(ind)+msc(10)
 447          nind=ind
 448        else
 449          DM(ind)=DM(ind)+isc(iabs(itype(nct)))
 450          DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
 451          nind=ind+1
 452         endif
 453         endif
 454
 455
 456         ind=1
 457         do i=nnt,nct
 458         mnum=molnum(i)
 459         if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
 460           .and.mnum.lt.4) then
 461           DU1(ind)=-isc(iabs(itype(i,1)))
 462           DU1(ind+1)=0.0d0
 463           ind=ind+2
 464         else
 465           DU1(ind)=mp/4-ip4
 466           ind=ind+1
 467         endif
 468         enddo
 469
 470        ind=1
 471        do i=nnt,nct-1
 472        mnum=molnum(i)
 473 !        if (iabs(itype(i,1)).eq.ntyp1) cycle
 474         write (iout,*) "i",i," itype",itype(i,1),ntyp1
 475         if (iabs(itype(i,1)).ne.10 .and. &
 476           iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
 477          DU2(ind)=mp/4-ip4
 478          DU2(ind+1)=0.0d0
 479          ind=ind+2
 480         else
 481          DU2(ind)=0.0d0
 482          DU2(ind+1)=0.0d0
 483          ind=ind+1
 484         endif
 485        enddo
 486 #else
 487       dimen=0
 488       dimen1=0
 489       dimenp=0
 490       do ichain=1,nchain
 491         dimen=dimen+chain_length(ichain)
 492         dimen1=dimen1+2*chain_length(ichain)-1
 493         dimenp=dimenp+chain_length(ichain)-1
 494       enddo
 495       write (iout,*) "Number of Calphas",dimen
 496       write (iout,*) "Number of sidechains",nside
 497       write (iout,*) "Number of peptide groups",dimenp
 498       dimen=dimen+nside ! number of centers
 499       dimen3=3*dimen  ! degrees of freedom
 500       write (iout,*) "Number of centers",dimen
 501       write (iout,*) "Degrees of freedom:",dimen3
 502 !      ip4=ip/4
 503       ind=1
 504       do ichain=1,nchain
 505         iposd_chain(ichain)=ind
 506         innt=chain_border(1,ichain)
 507         mnum=molnum(innt)
 508         inct=chain_border(2,ichain)
 509         if (mnum.eq.5) mp(mnum)=0.0
 510         if (mnum.eq.5) ip(mnum)=0.0
 511 !        if (mnum.eq.5) mp(mnum)=msc(itype(innt,mnum),mnum)
 512         DM(ind)=mp(mnum)/4+ip(mnum)/4
 513         if (iabs(itype(innt,1)).eq.10.or.molnum(innt).gt.2) then
 514            DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum)
 515           ind=ind+1
 516           nind=1
 517         else
 518           DM(ind)=DM(ind)+isc(iabs(itype(innt,mnum)),mnum)
 519           DM(ind+1)=msc(iabs(itype(innt,mnum)),mnum)+isc(iabs(itype(innt,mnum)),mnum)
 520           ind=ind+2
 521           nind=2
 522         endif
 523         write (iout,*) "ind",ind," nind",nind
 524         do i=innt+1,inct-1
 525          mnum=molnum(i)
 526 !        if (iabs(itype(i)).eq.ntyp1) cycle
 527 !          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 528 !          if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
 529         if (mnum.eq.5) mp(mnum)=0.0
 530         if (mnum.eq.5) ip(mnum)=0.0
 531 !        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 532           DM(ind)=2*ip(mnum)/4+mp(mnum)/2
 533           if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2) then
 534 !            if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)&
 535 !             DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum)
 536              DM(ind)=DM(ind)+msc(itype(i,mnum),mnum)
 537             ind=ind+1
 538             nind=nind+1
 539           else
 540             DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum)
 541             DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum)
 542             ind=ind+2
 543             nind=nind+2
 544           endif
 545           write (iout,*) "i",i," ind",ind," nind",nind
 546         enddo
 547         if (inct.gt.innt) then
 548 !         DM(ind)=ip4+mp(molnum(inct))/4
 549           mnum=molnum(inct)
 550         if (mnum.eq.5) mp(mnum)=0.0
 551 !        if (mnum.eq.5) mp(mnum)=msc(itype(inct,molnum(inct)),molnum(inct))
 552
 553           DM(ind)=mp(mnum)/4+ip(mnum)/4
 554           if (iabs(itype(inct,mnum)).eq.10.or.molnum(inct).gt.2) then
 555               DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct))
 556             ind=ind+1
 557             nind=nind+1
 558           else
 559             mnum=molnum(inct)
 560             DM(ind)=DM(ind)+isc(iabs(itype(inct,mnum)),mnum)
 561             DM(ind+1)=msc(iabs(itype(inct,mnum)),mnum)+isc(iabs(itype(inct,mnum)),mnum)
 562             ind=ind+2
 563             nind=nind+2
 564           endif
 565         endif
 566         write (iout,*) "ind",ind," nind",nind
 567         dimen_chain(ichain)=nind
 568         enddo
 569
 570       do ichain=1,nchain
 571         ind=iposd_chain(ichain)
 572         innt=chain_border(1,ichain)
 573         inct=chain_border(2,ichain)
 574         do i=innt,inct
 575          mnum=molnum(i)
 576           if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
 577             DU1(ind)=-isc(iabs(itype(i,mnum)),mnum)
 578             DU1(ind+1)=0.0d0
 579             ind=ind+2
 580           else
 581 !           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 582 !            if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
 583              if (mnum.eq.5) mp(mnum)=0.0
 584             DU1(ind)=mp(mnum)/4-ip(mnum)/4
 585             ind=ind+1
 586           endif
 587         enddo
 588       enddo
 589
 590       do ichain=1,nchain
 591         ind=iposd_chain(ichain)
 592         innt=chain_border(1,ichain)
 593         inct=chain_border(2,ichain)
 594         do i=innt,inct-1
 595          mnum=molnum(i)
 596 !       if (iabs(itype(i)).eq.ntyp1) cycle
 597 !c          write (iout,*) "i",i," itype",itype(i),ntyp1
 598           if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
 599 !            if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 600         if (mnum.eq.5) mp(mnum)=0.0
 601             DU2(ind)=mp(mnum)/4-ip(mnum)/4
 602             DU2(ind+1)=0.0d0
 603             ind=ind+2
 604           else
 605             DU2(ind)=0.0d0
 606             DU2(ind+1)=0.0d0
 607             ind=ind+1
 608           endif
 609         enddo
 610       enddo
 611       DMorig=DM
 612       DU1orig=DU1
 613       DU2orig=DU2
 614       gmatout=.true.
 615       if (gmatout) then
 616       write (iout,*)"The upper part of the five-diagonal inertia matrix"
 617       endif
 618       do ichain=1,nchain
 619         if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
 620         n=dimen_chain(ichain)
 621         innt=iposd_chain(ichain)
 622         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 623         if (gmatout) then
 624         do i=innt,inct
 625           if (i.lt.inct-1) then
 626             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 627           else if (i.eq.inct-1) then
 628             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 629           else
 630             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 631           endif
 632         enddo
 633         endif
 634         call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),&
 635         DU2(innt:inct-1), MARK)
 636
 637         if (mark.eq.-1) then
 638           write(iout,*)&
 639         "ERROR: the inertia matrix is not positive definite for chain",&
 640          ichain
 641 #ifdef MPI
 642          call MPI_Finalize(ierr)
 643 #endif
 644          stop
 645         else if (mark.eq.0) then
 646           write (iout,*)&
 647           "ERROR: the inertia matrix is singular for chain",ichain
 648 #ifdef MPI
 649           call MPI_Finalize(ierr)
 650 #endif
 651         else if (mark.eq.1) then
 652           if (gmatout) then
 653           write (iout,*) "The transformed five-diagonal inertia matrix"
 654           write (iout,'(a,i5)') 'Chain',ichain
 655           do i=innt,inct
 656             if (i.lt.inct-1) then
 657               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 658             else if (i.eq.inct-1) then
 659               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 660             else
 661               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 662             endif
 663           enddo
 664           endif
 665         endif
 666       enddo
 667 ! Diagonalization of the pentadiagonal matrix
 668 #ifdef TIMING
 669       time00=MPI_Wtime()
 670 #endif
 671
 672
 673 #endif
 674 !CRYSTBOND
 675 #else
 676       dimen=(nct-nnt+1)+nside
 677       dimen1=(nct-nnt)+(nct-nnt+1)
 678       dimen3=dimen*3
 679       write (iout,*) "nnt",nnt," nct",nct," nside",nside
 680 ! Old Gmatrix
 681 #ifdef MPI
 682       if (nfgtasks.gt.1) then
 683       time00=MPI_Wtime()
 684       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 685       time_Bcast=time_Bcast+MPI_Wtime()-time00
 686       call int_bounds(dimen,igmult_start,igmult_end)
 687       igmult_start=igmult_start-1
 688       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
 689           ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 690       my_ng_count=igmult_end-igmult_start
 691       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
 692           MPI_INTEGER,FG_COMM,IERROR)
 693       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
 694        ' absolute rank',myrank,' igmult_start',igmult_start,&
 695        ' igmult_end',igmult_end,' count',my_ng_count
 696       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 697       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 698       call flush(iout)
 699       else
 700 #endif
 701       igmult_start=1
 702       igmult_end=dimen
 703       my_ng_count=dimen
 704 #ifdef MPI
 705       endif
 706 #endif
 707 !      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 708 !  Zeroing out A and fricmat
 709       do i=1,dimen
 710         do j=1,dimen
 711           A(i,j)=0.0D0
 712         enddo
 713       enddo
 714 !  Diagonal elements of the dC part of A and the respective friction coefficients
 715       ind=1
 716       ind1=0
 717 !      print *,"TUTUTUT?!",nnt,nct-1
 718       do i=nnt,nct-1
 719         mnum=molnum(i)
 720         ind=ind+1
 721         ind1=ind1+1
 722         coeff=0.25d0*IP(mnum)
 723         print *,"TU",i,itype(i,mnum),mnum
 724         if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
 725         print *,"TU2",i,coeff,mp(mnum)
 726         massvec(ind1)=mp(mnum)
 727         Gmat(ind,ind)=coeff
 728 !        print *,"i",mp(mnum)
 729         A(ind1,ind)=0.5d0
 730       enddo
 731
 732 !  Off-diagonal elements of the dC part of A
 733       k=3
 734       do i=1,nct-nnt
 735         do j=1,i
 736           A(i,j)=1.0d0
 737         enddo
 738       enddo
 739 !  Diagonal elements of the dX part of A and the respective friction coefficients
 740       m=nct-nnt
 741       m1=nct-nnt+1
 742       ind=0
 743       ind1=0
 744       do i=1,2
 745       msc(ntyp1_molec(i),i)=1.0d0
 746       enddo
 747       msc(ntyp1_molec(4),4)=1.0d0
 748 !      print *,"TU3",ntyp1_molec(4)-1
 749        msc(ntyp1_molec(4)-1,4)=1.0d0
 750       do i=nnt,nct
 751         ind=ind+1
 752         ii = ind+m
 753         mnum=molnum(i)
 754         iti=itype(i,mnum)
 755 !        if (mnum.eq.5) then
 756 !        mscab=0.0
 757 !        else
 758         mscab=msc(iabs(iti),mnum)
 759 !        endif
 760         massvec(ii)=mscab
 761
 762         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
 763           ind1=ind1+1
 764           ii1= ind1+m1
 765           A(ii,ii1)=1.0d0
 766           Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
 767         endif
 768       enddo
 769 !  Off-diagonal elements of the dX part of A
 770       ind=0
 771       k=nct-nnt
 772       do i=nnt,nct
 773         mnum=molnum(i)
 774         iti=itype(i,mnum)
 775         ind=ind+1
 776         do j=nnt,i
 777           ii = ind
 778           jj = j-nnt+1
 779           A(k+ii,jj)=1.0d0
 780         enddo
 781       enddo
 782       if (lprn) then
 783         write (iout,*)
 784         write (iout,*) "Vector massvec"
 785         do i=1,dimen1
 786           write (iout,*) i,massvec(i)
 787         enddo
 788         write (iout,'(//a)') "A"
 789         call matout(dimen,dimen1,nres2,nres2,A)
 790       endif
 791
 792 ! Calculate the G matrix (store in Gmat)
 793       do k=1,dimen
 794        do i=1,dimen
 795          dtdi=0.0d0
 796          do j=1,dimen1
 797            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 798          enddo
 799          Gmat(k,i)=Gmat(k,i)+dtdi
 800        enddo
 801       enddo
 802
 803       if (lprn) then
 804         write (iout,'(//a)') "Gmat"
 805         call matout(dimen,dimen,nres2,nres2,Gmat)
 806       endif
 807       do i=1,dimen
 808         do j=1,dimen
 809           Ginv(i,j)=0.0d0
 810           Gcopy(i,j)=Gmat(i,j)
 811         enddo
 812         Ginv(i,i)=1.0d0
 813       enddo
 814 ! Invert the G matrix
 815       call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
 816       if (lprn) then
 817         write (iout,'(//a)') "Ginv"
 818         call matout(dimen,dimen,nres2,nres2,Ginv)
 819       endif
 820
 821       Bmat=0.0d0
 822       Bmat(1,1)=1.0d0
 823       j=2
 824       do i=nnt,nct
 825         mnum=molnum(i)
 826         if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
 827           .or.(mnum.ge.4)) then
 828           Bmat(i-nnt+2,j-1)=-1
 829           Bmat(i-nnt+2,j)=1
 830           j=j+1
 831         else
 832           if (i.lt.nct) then
 833             Bmat(i-nnt+2,j-1)=-1
 834             Bmat(i-nnt+2,j+1)=1
 835             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 836             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 837             j=j+2
 838           else
 839             Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
 840             Bmat(i-nnt+1+(nct-nnt)+1,j)=1
 841           endif
 842         endif
 843       enddo
 844       write (iout,*) "j",j," dimen",dimen
 845       if (lprn) then
 846         write (iout,'(//a)') "Bmat"
 847         call matout(dimen,dimen,nres2,nres2,Bmat)
 848       endif
 849       Gcopy=0.0d0
 850       write(iout,*) "before Gcopy",dimen,nres*2
 851 #ifdef TESTFIVEDIAG
 852       do i=1,dimen
 853         do j=1,dimen
 854 !            write (iout,*) "myij",i,j
 855             do k=1,dimen
 856                do l=1,dimen
 857 !                 write(iout,*) "i,j,k,l",i,j,k,l
 858                  Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
 859                enddo
 860             enddo
 861         enddo
 862       enddo
 863 #endif
 864       if (lprn) then
 865         write (iout,'(//a)') "Gmat-transformed"
 866         call matout(dimen,dimen,nres2,nres2,Gcopy)
 867       endif
 868 #ifdef MPI
 869       if (nfgtasks.gt.1) then
 870         myginv_ng_count=nres2*my_ng_count
 871         call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
 872           nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 873         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
 874           nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 875         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 876         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 877         call flush(iout)
 878 !        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 879 !     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 880 !     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 881 !        call MPI_Barrier(FG_COMM,IERR)
 882         time00=MPI_Wtime()
 883         call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
 884           nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
 885           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 886 #ifdef TIMING
 887         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 888 #endif
 889         do i=1,dimen
 890           do j=1,2*my_ng_count
 891             ginv(j,i)=gcopy(i,j)
 892           enddo
 893         enddo
 894 !        write (iout,*) "Master's chunk of ginv"
 895 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 896       endif
 897 #endif
 898       if (osob) then
 899         write (iout,*) "The G matrix is singular."
 900         write (iout,'(//a)') "Gmat-transformed"
 901         call matout(dimen,dimen,nres2,nres2,Gcopy)
 902         stop
 903       endif
 904 ! Compute G**(-1/2) and G**(1/2)
 905       ind=0
 906       do i=1,dimen
 907         do j=1,i
 908           ind=ind+1
 909           Ghalf(ind)=Gmat(i,j)
 910         enddo
 911       enddo
 912       call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
 913         ierr,iwork)
 914       if (lprn) then
 915         write (iout,'(//a)') &
 916          "Eigenvectors and eigenvalues of the G matrix"
 917         call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
 918       endif
 919       do i=1,dimen
 920         sqreig(i)=dsqrt(Geigen(i))
 921       enddo
 922       do i=1,dimen
 923         do j=1,dimen
 924           Gsqrp(i,j)=0.0d0
 925           Gsqrm(i,j)=0.0d0
 926           Gcopy(i,j)=0.0d0
 927           do k=1,dimen
 928             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 929             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 930             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 931           enddo
 932         enddo
 933       enddo
 934       if (lprn) then
 935         write (iout,*) "Comparison of original and restored G"
 936         do i=1,dimen
 937           do j=1,dimen
 938             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
 939               Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 940           enddo
 941         enddo
 942       endif
 943       if (allocated(Gcopy)) deallocate(Gcopy)
 944 #endif
 945 !write(iout,*) "end setup_MD_matr"
 946       return
 947       end subroutine setup_MD_matrices
 948 !-----------------------------------------------------------------------------
 949       subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
 950
 951 !      implicit real*8 (a-h,o-z)
 952 !      include 'DIMENSIONS'
 953 !      include 'COMMON.IOUNITS'
 954       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 955       real(kind=8) :: A(LM2,LM3),B(LM2)
 956       KA=1
 957       KC=6
 958     1 KB=MIN0(KC,NC)
 959       WRITE(IOUT,600) (I,I=KA,KB)
 960       WRITE(IOUT,601) (B(I),I=KA,KB)
 961       WRITE(IOUT,602)
 962     2 N=0
 963       DO 3  I=1,NR
 964       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 965       N=N+1
 966       IF(N.LT.10) GO TO 3
 967       WRITE(IOUT,602)
 968       N=0
 969     3 CONTINUE
 970     4 IF (KB.EQ.NC) RETURN
 971       KA=KC+1
 972       KC=KC+6
 973       GO TO 1
 974   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 975   601 FORMAT (/5X,10(1PE11.4))
 976   602 FORMAT (2H  )
 977   603 FORMAT (I5,10F11.5)
 978   604 FORMAT (1H1)
 979       end subroutine EIGOUT
 980 !-----------------------------------------------------------------------------
 981       subroutine MATOUT(NC,NR,LM2,LM3,A)
 982
 983 !      implicit real*8 (a-h,o-z)
 984 !      include 'DIMENSIONS'
 985 !      include 'COMMON.IOUNITS'
 986       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 987       real(kind=8) :: A(LM2,LM3)
 988       KA=1
 989       KC=6
 990     1 KB=MIN0(KC,NC)
 991       WRITE(IOUT,600) (I,I=KA,KB)
 992       WRITE(IOUT,602)
 993     2 N=0
 994       DO 3  I=1,NR
 995       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 996       N=N+1
 997       IF(N.LT.10) GO TO 3
 998       WRITE(IOUT,602)
 999       N=0
1000     3 CONTINUE
1001     4 IF (KB.EQ.NC) RETURN
1002       KA=KC+1
1003       KC=KC+6
1004       GO TO 1
1005   600 FORMAT (//5x,9I11)
1006   602 FORMAT (2H  )
1007   603 FORMAT (I5,10F11.3)
1008   604 FORMAT (1H1)
1009       end subroutine MATOUT
1010 !-----------------------------------------------------------------------------
1011       subroutine MATOUT1(NC,NR,LM2,LM3,A)
1012
1013 !      implicit real*8 (a-h,o-z)
1014 !      include 'DIMENSIONS'
1015 !      include 'COMMON.IOUNITS'
1016       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
1017       real(kind=8) :: A(LM2,LM3)
1018       KA=1
1019       KC=21
1020     1 KB=MIN0(KC,NC)
1021       WRITE(IOUT,600) (I,I=KA,KB)
1022       WRITE(IOUT,602)
1023     2 N=0
1024       DO 3  I=1,NR
1025       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
1026       N=N+1
1027       IF(N.LT.3) GO TO 3
1028       WRITE(IOUT,602)
1029       N=0
1030     3 CONTINUE
1031     4 IF (KB.EQ.NC) RETURN
1032       KA=KC+1
1033       KC=KC+21
1034       GO TO 1
1035   600 FORMAT (//5x,7(3I5,2x))
1036   602 FORMAT (2H  )
1037   603 FORMAT (I5,7(3F5.1,2x))
1038   604 FORMAT (1H1)
1039       end subroutine MATOUT1
1040 !-----------------------------------------------------------------------------
1041       subroutine MATOUT2(NC,NR,LM2,LM3,A)
1042
1043 !      implicit real*8 (a-h,o-z)
1044 !      include 'DIMENSIONS'
1045 !      include 'COMMON.IOUNITS'
1046       integer :: I,J,KA,KC,KB,N
1047       integer :: LM2,LM3,NC,NR
1048       real(kind=8) :: A(LM2,LM3)
1049       KA=1
1050       KC=12
1051     1 KB=MIN0(KC,NC)
1052       WRITE(IOUT,600) (I,I=KA,KB)
1053       WRITE(IOUT,602)
1054     2 N=0
1055       DO 3  I=1,NR
1056       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
1057       N=N+1
1058       IF(N.LT.3) GO TO 3
1059       WRITE(IOUT,602)
1060       N=0
1061     3 CONTINUE
1062     4 IF (KB.EQ.NC) RETURN
1063       KA=KC+1
1064       KC=KC+12
1065       GO TO 1
1066   600 FORMAT (//5x,4(3I9,2x))
1067   602 FORMAT (2H  )
1068   603 FORMAT (I5,4(3F9.3,2x))
1069   604 FORMAT (1H1)
1070       end subroutine MATOUT2
1071 #ifdef FIVEDIAG
1072       subroutine fivediagmult(n,DM,DU1,DU2,x,y)
1073       integer n
1074       double precision DM(n),DU1(n),DU2(n),x(n),y(n)
1075       integer i
1076       y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
1077       y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
1078       do i=3,n-2
1079         y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i) &
1080            +DU1(i)*x(i+1)+DU2(i)*x(i+2)
1081       enddo
1082       y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1) &
1083       +DU1(n-1)*x(n)
1084       y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
1085       return
1086       end subroutine
1087 !c---------------------------------------------------------------------------
1088       subroutine fivediaginv_mult(ndim,forces,d_a_vec)
1089       use energy_data, only:nchain,chain_border,nct,nnt,molnum,&
1090       chain_border1,itype
1091       use geometry_data, only: nside
1092       integer ndim
1093
1094       real(kind=8),dimension(:),allocatable ::forces,d_a_vec
1095       real(kind=8),dimension(:),allocatable :: xsolv,rs
1096       real(kind=8),dimension(:,:),allocatable :: accel
1097       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev,mnum
1098       accel=0.0d0
1099 !#define DEBUG
1100       if (.not.allocated(forces)) allocate(forces(6*nres))
1101       if (.not.allocated(d_a_vec)) allocate(d_a_vec(6*nres))
1102       if (.not.allocated(xsolv)) allocate(xsolv(3*ndim))
1103       if (.not.allocated(rs)) allocate(rs(3*ndim))
1104       if (.not.allocated(accel)) allocate(accel(3,0:2*nres))
1105
1106 #ifdef DEBUG
1107       do i=1,6*nres
1108        write(iout,*) "received forces",i,forces(i)
1109       enddo
1110 #endif
1111       do j=1,3
1112 !Compute accelerations in Calpha and SC
1113         do ichain=1,nchain
1114           n=dimen_chain(ichain)
1115           iposc=iposd_chain(ichain)
1116           innt=chain_border(1,ichain)
1117           inct=chain_border(2,ichain)
1118           do i=iposc,iposc+n-1
1119 !            write(iout,*) "index",3*(i-1)+j,forces(3*(i-1)+j)
1120             rs(i-iposc+1)=forces(3*(i-1)+j)
1121           enddo
1122 #ifdef DEBUG
1123           write (iout,*) "j",j," chain",ichain,"n",n
1124           write (iout,*) "innt",innt,inct,iposc
1125           write (iout,*) "rs"
1126           do i=1,n
1127           write (iout,'(i5,f10.5)') i,rs(i)
1128           enddo
1129 #endif
1130           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
1131 #ifdef DEBUG
1132           write (iout,*) "xsolv"
1133           write (iout,'(f10.5)') (xsolv(i),i=1,n)
1134 #endif
1135           ind=1
1136           do i=innt,inct
1137             mnum=molnum(i)
1138             if (itype(i,1).eq.10.or.mnum.gt.2)then
1139               accel(j,i)=xsolv(ind)
1140               ind=ind+1
1141             else
1142               accel(j,i)=xsolv(ind)
1143               accel(j,i+nres)=xsolv(ind+1)
1144               ind=ind+2
1145             end if
1146           enddo
1147         enddo
1148       enddo
1149 !C Convert d_a to virtual-bon-vector basis
1150 #ifdef DEBUG
1151       write (iout,*) "accel in CA-SC basis"
1152       do i=1,nres
1153         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),&
1154            (accel(j,i+nres),j=1,3)
1155       enddo
1156       write (iout,*) "nnt",nnt
1157 #endif
1158       if (nnt.eq.1) then
1159         accel(:,0)=accel(:,1)
1160       endif
1161       do i=1,nres
1162       mnum=molnum(i)
1163         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
1164         .or.mnum.ge.3) then
1165           do j=1,3
1166             accel(j,i)=accel(j,i+1)-accel(j,i)
1167           enddo
1168         else
1169           do j=1,3
1170             accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
1171             accel(j,i)=accel(j,i+1)-accel(j,i)
1172           enddo
1173         end if
1174       enddo
1175       accel(:,nres)=0.0d0
1176       accel(:,nct)=0.0d0
1177       accel(:,2*nres)=0.0d0
1178       if (nnt.gt.1) then
1179         accel(:,0)=accel(:,1)
1180         accel(:,1)=0.0d0
1181       endif
1182       do ichain=2,nchain
1183       if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
1184         accel(:,chain_border1(1,ichain)-1)= &
1185          accel(:,chain_border1(1,ichain))
1186         accel(:,chain_border1(1,ichain))=0.0d0
1187       enddo
1188       do ichain=2,nchain
1189       if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
1190         accel(:,chain_border1(1,ichain)-1)= &
1191        accel(:,chain_border1(1,ichain)-1) &
1192         +accel(:,chain_border(2,ichain-1))
1193         accel(:,chain_border(2,ichain-1))=0.0d0
1194       enddo
1195 #ifdef DEBUG
1196       write (iout,*) "accel in fivediaginv_mult: 1"
1197       do i=0,2*nres
1198         write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
1199       enddo
1200 #endif
1201       do j=1,3
1202         d_a_vec(j)=accel(j,0)
1203       enddo
1204       ind=3
1205       do i=nnt,nct-1
1206         do j=1,3
1207           d_a_vec(ind+j)=accel(j,i)
1208         enddo
1209         ind=ind+3
1210       enddo
1211       do i=nnt,nct
1212         mnum=molnum(i)
1213         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1214           .and.mnum.le.2) then
1215           do j=1,3
1216             d_a_vec(ind+j)=accel(j,i+nres)
1217           enddo
1218           ind=ind+3
1219         endif
1220       enddo
1221 #ifdef DEBUG
1222       write (iout,*) "d_a_vec"
1223       write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
1224 #endif
1225       return
1226       end subroutine
1227 #undef DEBUG
1228
1229 #else
1230 !-----------------------------------------------------------------------------
1231       subroutine ginv_mult(z,d_a_tmp)
1232
1233       use geometry_data, only: nres
1234       use control_data
1235       use MPI_data
1236 !      implicit real*8 (a-h,o-z)
1237 !      include 'DIMENSIONS'
1238 #ifdef MPI
1239       include 'mpif.h'
1240       integer :: ierr,ierror
1241 #endif
1242 !      include 'COMMON.SETUP'
1243 !      include 'COMMON.TIME1'
1244 !      include 'COMMON.MD'
1245       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1246       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1247       real(kind=8) :: time00,time01
1248       integer :: i,j,k,ind
1249 #ifdef MPI
1250       if (nfgtasks.gt.1) then
1251         if (fg_rank.eq.0) then
1252 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1253           time00=MPI_Wtime()
1254           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
1255           time_Bcast=time_Bcast+MPI_Wtime()-time00
1256 !          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
1257         endif
1258 !        write (2,*) "time00",time00
1259 !        write (2,*) "Before Scatterv"
1260 !        call flush(2)
1261 !        write (2,*) "Whole z (for FG master)"
1262 !        do i=1,dimen
1263 !          write (2,*) i,z(i)
1264 !        enddo
1265 !        call MPI_Barrier(FG_COMM,IERROR)
1266         time00=MPI_Wtime()
1267 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
1268         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1269           MPI_DOUBLE_PRECISION,&
1270           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1271 !        write (2,*) "My chunk of z"
1272 !        do i=1,3*my_ng_count
1273 !          write (2,*) i,z(i)
1274 !        enddo
1275 !        write (2,*) "After SCATTERV"
1276 !        call flush(2)
1277 !        write (2,*) "MPI_Wtime",MPI_Wtime()
1278         time_scatter=time_scatter+MPI_Wtime()-time00
1279 #ifdef TIMING
1280         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1281 #endif
1282 !        write (2,*) "time_scatter",time_scatter
1283 !        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1284 !     &    my_ng_count
1285 !        call flush(2)
1286         time01=MPI_Wtime()
1287         do k=0,2
1288           do i=1,dimen
1289             ind=(i-1)*3+k+1
1290             temp(ind)=0.0d0
1291             do j=1,my_ng_count
1292 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1293 !     &         Ginv(i,j),z((j-1)*3+k+1),
1294 !     &          Ginv(i,j)*z((j-1)*3+k+1)
1295 !              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1296               temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
1297             enddo
1298           enddo
1299         enddo
1300         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1301 !        write (2,*) "Before REDUCE"
1302 !        call flush(2)
1303 !        write (2,*) "z before reduce"
1304 !        do i=1,dimen
1305 !          write (2,*) i,temp(i)
1306 !        enddo
1307         time00=MPI_Wtime()
1308         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1309             MPI_SUM,king,FG_COMM,IERR)
1310         time_reduce=time_reduce+MPI_Wtime()-time00
1311 !        write (2,*) "After REDUCE"
1312 !        call flush(2)
1313       else
1314 #endif
1315 #ifdef TIMING
1316         time01=MPI_Wtime()
1317 #endif
1318         do k=0,2
1319           do i=1,dimen
1320             ind=(i-1)*3+k+1
1321             d_a_tmp(ind)=0.0d0
1322             do j=1,dimen
1323 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1324 !              call flush(2)
1325 !     &         Ginv(i,j),z((j-1)*3+k+1),
1326 !     &          Ginv(i,j)*z((j-1)*3+k+1)
1327               d_a_tmp(ind)=d_a_tmp(ind) &
1328                                +Ginv(j,i)*z((j-1)*3+k+1)
1329 !              d_a_tmp(ind)=d_a_tmp(ind)
1330 !     &                         +Ginv(i,j)*z((j-1)*3+k+1)
1331             enddo
1332           enddo
1333         enddo
1334 #ifdef TIMING
1335         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1336 #endif
1337 #ifdef MPI
1338       endif
1339 #endif
1340       return
1341       end subroutine ginv_mult
1342 !-----------------------------------------------------------------------------
1343 #ifdef GINV_MULT
1344       subroutine ginv_mult_test(z,d_a_tmp)
1345
1346 !      include 'DIMENSIONS'
1347 !el      integer :: dimen
1348 !      include 'COMMON.MD'
1349       real(kind=8),dimension(dimen) :: z,d_a_tmp
1350       real(kind=8),dimension(dimen/3) :: ztmp,dtmp
1351       integer :: i,j,k,ind
1352 !      do i=1,dimen
1353 !        d_a_tmp(i)=0.0d0
1354 !        do j=1,dimen
1355 !          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1356 !        enddo
1357 !      enddo
1358 !
1359 !      return
1360
1361 !ibm* unroll(3)
1362       do k=0,2
1363        do j=1,dimen/3
1364         ztmp(j)=z((j-1)*3+k+1)
1365        enddo
1366
1367        call alignx(16,ztmp(1))
1368        call alignx(16,dtmp(1))
1369        call alignx(16,Ginv(1,1))
1370
1371        do i=1,dimen/3
1372         dtmp(i)=0.0d0
1373         do j=1,dimen/3
1374            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1375         enddo
1376        enddo
1377        do i=1,dimen/3
1378         ind=(i-1)*3+k+1
1379         d_a_tmp(ind)=dtmp(i)
1380        enddo
1381       enddo
1382       return
1383       end subroutine ginv_mult_test
1384 #endif
1385 !-----------------------------------------------------------------------------
1386       subroutine fricmat_mult(z,d_a_tmp)
1387
1388       use geometry_data, only: nres
1389       use control_data
1390       use MPI_data
1391 !      include 'DIMENSIONS'
1392 #ifdef MPI
1393       include 'mpif.h'
1394       integer :: IERROR,ierr
1395 #endif
1396 !      include 'COMMON.MD'
1397 !      include 'COMMON.IOUNITS'
1398 !      include 'COMMON.SETUP'
1399 !      include 'COMMON.TIME1'
1400 !#ifndef LANG0
1401 !      include 'COMMON.LANGEVIN'
1402 !#else
1403 !      include 'COMMON.LANGEVIN.lang0'
1404 !#endif
1405       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1406       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
1407       real(kind=8) :: time00,time01
1408       integer :: i,j,k,ind,nres2
1409       nres2=2*nres
1410 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1411
1412 #ifdef MPI
1413       if (nfgtasks.gt.1) then
1414         if (fg_rank.eq.0) then
1415 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1416           time00=MPI_Wtime()
1417           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1418           time_Bcast=time_Bcast+MPI_Wtime()-time00
1419 !          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1420         endif
1421 !        call MPI_Barrier(FG_COMM,IERROR)
1422         time00=MPI_Wtime()
1423         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1424           MPI_DOUBLE_PRECISION,&
1425           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1426 !        write (2,*) "My chunk of z"
1427 !        do i=1,3*my_ng_count
1428 !          write (2,*) i,z(i)
1429 !        enddo
1430         time_scatter=time_scatter+MPI_Wtime()-time00
1431 #ifdef TIMING
1432         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1433 #endif
1434         time01=MPI_Wtime()
1435         do k=0,2
1436           do i=1,dimen
1437             ind=(i-1)*3+k+1
1438             temp(ind)=0.0d0
1439             do j=1,my_ng_count
1440               temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1441             enddo
1442           enddo
1443         enddo
1444         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1445 !        write (2,*) "Before REDUCE"
1446 !        write (2,*) "d_a_tmp before reduce"
1447 !        do i=1,dimen3
1448 !          write (2,*) i,temp(i)
1449 !        enddo
1450 !        call flush(2)
1451         time00=MPI_Wtime()
1452         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1453             MPI_SUM,king,FG_COMM,IERR)
1454         time_reduce=time_reduce+MPI_Wtime()-time00
1455 !        write (2,*) "After REDUCE"
1456 !        call flush(2)
1457       else
1458 #endif
1459 #ifdef TIMING
1460         time01=MPI_Wtime()
1461 #endif
1462         do k=0,2
1463          do i=1,dimen
1464           ind=(i-1)*3+k+1
1465           d_a_tmp(ind)=0.0d0
1466           do j=1,dimen
1467              d_a_tmp(ind)=d_a_tmp(ind) &
1468                                  -fricmat(j,i)*z((j-1)*3+k+1)
1469           enddo
1470          enddo
1471         enddo
1472 #ifdef TIMING
1473         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1474 #endif
1475 #ifdef MPI
1476       endif
1477 #endif
1478 !      write (iout,*) "Vector d_a"
1479 !      do i=1,dimen3
1480 !        write (2,*) i,d_a_tmp(i)
1481 !      enddo
1482       return
1483       end subroutine fricmat_mult
1484 !-----------------------------------------------------------------------------
1485 #endif
1486 !-----------------------------------------------------------------------------
1487       end module REMD