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