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