source/unres/src_MD-M/lagrangian_lesyng.F

   1        subroutine lagrangian
   2 c-------------------------------------------------------------------------
   3 c  This subroutine contains the total lagrangain from which the accelerations
   4 c  are obtained.  For numerical gradient checking, the derivetive of the
   5 c  lagrangian in the velocities and coordinates are calculated seperately
   6 c-------------------------------------------------------------------------
   7        implicit real*8 (a-h,o-z)
   8        include 'DIMENSIONS'
   9 #ifdef MPI
  10        include 'mpif.h'
  11 #endif
  12        include 'COMMON.VAR'
  13        include 'COMMON.CHAIN'
  14        include 'COMMON.DERIV'
  15        include 'COMMON.GEO'
  16        include 'COMMON.LOCAL'
  17        include 'COMMON.INTERACT'
  18        include 'COMMON.MD'
  19        include 'COMMON.IOUNITS'
  20        include 'COMMON.CONTROL'
  21        include 'COMMON.MUCA'
  22        include 'COMMON.TIME1'
  23
  24        integer i,j,ind
  25        double precision zapas(MAXRES6),muca_factor
  26        logical lprn /.false./
  27        common /cipiszcze/ itime
  28
  29 #ifdef TIMING
  30        time00=MPI_Wtime()
  31 #endif
  32        do j=1,3
  33          zapas(j)=-gcart(j,0)
  34        enddo
  35       ind=3
  36       if (lprn) then
  37         write (iout,*) "Potential forces backbone"
  38       endif
  39       do i=nnt,nct-1
  40         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
  41      &    i,(-gcart(j,i),j=1,3)
  42         do j=1,3
  43           ind=ind+1
  44           zapas(ind)=-gcart(j,i)
  45         enddo
  46       enddo
  47       if (lprn) write (iout,*) "Potential forces sidechain"
  48       do i=nnt,nct
  49         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
  50           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
  51      &       i,(-gcart(j,i),j=1,3)
  52           do j=1,3
  53             ind=ind+1
  54             zapas(ind)=-gxcart(j,i)
  55           enddo
  56         endif
  57       enddo
  58
  59       call ginv_mult(zapas,d_a_work)
  60
  61       do j=1,3
  62         d_a(j,0)=d_a_work(j)
  63       enddo
  64       ind=3
  65       do i=nnt,nct-1
  66         do j=1,3
  67           ind=ind+1
  68           d_a(j,i)=d_a_work(ind)
  69         enddo
  70       enddo
  71       do i=nnt,nct
  72         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
  73           do j=1,3
  74             ind=ind+1
  75             d_a(j,i+nres)=d_a_work(ind)
  76           enddo
  77         endif
  78       enddo
  79
  80       if(lmuca) then
  81        imtime=imtime+1
  82        if(mucadyn.gt.0) call muca_update(potE)
  83        factor=muca_factor(potE)*t_bath*Rb
  84
  85 cd       print *,'lmuca ',factor,potE
  86        do j=1,3
  87           d_a(j,0)=d_a(j,0)*factor
  88        enddo
  89        do i=nnt,nct-1
  90          do j=1,3
  91           d_a(j,i)=d_a(j,i)*factor
  92          enddo
  93        enddo
  94        do i=nnt,nct
  95          do j=1,3
  96           d_a(j,i+nres)=d_a(j,i+nres)*factor
  97          enddo
  98        enddo
  99
 100       endif
 101
 102       if (lprn) then
 103         write(iout,*) 'acceleration 3D'
 104         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 105         do i=nnt,nct-1
 106          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 107         enddo
 108         do i=nnt,nct
 109          write (iout,'(i3,3f10.5,3x,3f10.5)')
 110      &     i+nres,(d_a(j,i+nres),j=1,3)
 111         enddo
 112       endif
 113 #ifdef TIMING
 114       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 115 #endif
 116       return
 117       end
 118 c------------------------------------------------------------------
 119       subroutine setup_MD_matrices
 120       implicit real*8 (a-h,o-z)
 121       include 'DIMENSIONS'
 122 #ifdef MPI
 123       include 'mpif.h'
 124       integer ierror
 125 #endif
 126       include 'COMMON.SETUP'
 127       include 'COMMON.VAR'
 128       include 'COMMON.CHAIN'
 129       include 'COMMON.DERIV'
 130       include 'COMMON.GEO'
 131       include 'COMMON.LOCAL'
 132       include 'COMMON.INTERACT'
 133       include 'COMMON.MD'
 134 #ifndef LANG0
 135       include 'COMMON.LANGEVIN'
 136 #else
 137       include 'COMMON.LANGEVIN.lang0'
 138 #endif
 139       include 'COMMON.IOUNITS'
 140       include 'COMMON.TIME1'
 141       integer i,j
 142       logical lprn /.false./
 143       logical osob
 144       double precision dtdi,massvec(maxres2),Gcopy(maxres2,maxres2),
 145      &  Ghalf(mmaxres2),sqreig(maxres2)
 146       double precision work(8*maxres6)
 147       integer iwork(maxres6)
 148       common /przechowalnia/ Gcopy,Ghalf
 149 c
 150 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 151 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 152 c
 153 c Determine the number of degrees of freedom (dimen) and the number of
 154 c sites (dimen1)
 155       dimen=(nct-nnt+1)+nside
 156       dimen1=(nct-nnt)+(nct-nnt+1)
 157       dimen3=dimen*3
 158 #ifdef MPI
 159       if (nfgtasks.gt.1) then
 160       time00=MPI_Wtime()
 161       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 162       time_Bcast=time_Bcast+MPI_Wtime()-time00
 163       call int_bounds(dimen,igmult_start,igmult_end)
 164       igmult_start=igmult_start-1
 165       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 166      &    ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 167       my_ng_count=igmult_end-igmult_start
 168       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 169      &    MPI_INTEGER,FG_COMM,IERROR)
 170       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
 171      & ' absolute rank',myrank,' igmult_start',igmult_start,
 172      & ' igmult_end',igmult_end,' count',my_ng_count
 173       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 174       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 175       call flush(iout)
 176       else
 177 #endif
 178       igmult_start=1
 179       igmult_end=dimen
 180       my_ng_count=dimen
 181 #ifdef MPI
 182       endif
 183 #endif
 184 c      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 185 c  Zeroing out A and fricmat
 186       do i=1,dimen
 187         do j=1,dimen
 188           A(i,j)=0.0D0
 189         enddo
 190       enddo
 191 c  Diagonal elements of the dC part of A and the respective friction coefficients
 192       ind=1
 193       ind1=0
 194       do i=nnt,nct-1
 195         ind=ind+1
 196         ind1=ind1+1
 197         coeff=0.25d0*IP
 198         massvec(ind1)=mp
 199         Gmat(ind,ind)=coeff
 200         A(ind1,ind)=0.5d0
 201       enddo
 202
 203 c  Off-diagonal elements of the dC part of A
 204       k=3
 205       do i=1,nct-nnt
 206         do j=1,i
 207           A(i,j)=1.0d0
 208         enddo
 209       enddo
 210 c  Diagonal elements of the dX part of A and the respective friction coefficients
 211       m=nct-nnt
 212       m1=nct-nnt+1
 213       ind=0
 214       ind1=0
 215       msc(ntyp1)=1.0d0
 216       do i=nnt,nct
 217         ind=ind+1
 218         ii = ind+m
 219         iti=itype(i)
 220         massvec(ii)=msc(iabs(iti))
 221         if (iti.ne.10 .and. iti.ne.ntyp1) then
 222           ind1=ind1+1
 223           ii1= ind1+m1
 224           A(ii,ii1)=1.0d0
 225           Gmat(ii1,ii1)=ISC(iabs(iti))
 226         endif
 227       enddo
 228 c  Off-diagonal elements of the dX part of A
 229       ind=0
 230       k=nct-nnt
 231       do i=nnt,nct
 232         iti=itype(i)
 233         ind=ind+1
 234         do j=nnt,i
 235           ii = ind
 236           jj = j-nnt+1
 237           A(k+ii,jj)=1.0d0
 238         enddo
 239       enddo
 240       if (lprn) then
 241         write (iout,*)
 242         write (iout,*) "Vector massvec"
 243         do i=1,dimen1
 244           write (iout,*) i,massvec(i)
 245         enddo
 246         write (iout,'(//a)') "A"
 247         call matout(dimen,dimen1,maxres2,maxres2,A)
 248       endif
 249
 250 c Calculate the G matrix (store in Gmat)
 251       do k=1,dimen
 252        do i=1,dimen
 253          dtdi=0.0d0
 254          do j=1,dimen1
 255            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 256          enddo
 257          Gmat(k,i)=Gmat(k,i)+dtdi
 258        enddo
 259       enddo
 260
 261       if (lprn) then
 262         write (iout,'(//a)') "Gmat"
 263         call matout(dimen,dimen,maxres2,maxres2,Gmat)
 264       endif
 265       do i=1,dimen
 266         do j=1,dimen
 267           Ginv(i,j)=0.0d0
 268           Gcopy(i,j)=Gmat(i,j)
 269         enddo
 270         Ginv(i,i)=1.0d0
 271       enddo
 272 c Invert the G matrix
 273       call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
 274       if (lprn) then
 275         write (iout,'(//a)') "Ginv"
 276         call matout(dimen,dimen,maxres2,maxres2,Ginv)
 277       endif
 278 #ifdef MPI
 279       if (nfgtasks.gt.1) then
 280         myginv_ng_count=maxres2*my_ng_count
 281         call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 282      &    nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 283         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 284      &    nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 285         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 286         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 287         call flush(iout)
 288 c        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 289 c     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 290 c     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 291 c        call MPI_Barrier(FG_COMM,IERR)
 292         time00=MPI_Wtime()
 293         call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 294      &    nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 295      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 296 #ifdef TIMING
 297         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 298 #endif
 299         do i=1,dimen
 300           do j=1,2*my_ng_count
 301             ginv(j,i)=gcopy(i,j)
 302           enddo
 303         enddo
 304 c        write (iout,*) "Master's chunk of ginv"
 305 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 306       endif
 307 #endif
 308       if (osob) then
 309         write (iout,*) "The G matrix is singular."
 310         stop
 311       endif
 312 c Compute G**(-1/2) and G**(1/2)
 313       ind=0
 314       do i=1,dimen
 315         do j=1,i
 316           ind=ind+1
 317           Ghalf(ind)=Gmat(i,j)
 318         enddo
 319       enddo
 320       call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
 321      &  ierr,iwork)
 322       if (lprn) then
 323         write (iout,'(//a)')
 324      &   "Eigenvectors and eigenvalues of the G matrix"
 325         call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
 326       endif
 327       do i=1,dimen
 328         sqreig(i)=dsqrt(Geigen(i))
 329       enddo
 330       do i=1,dimen
 331         do j=1,dimen
 332           Gsqrp(i,j)=0.0d0
 333           Gsqrm(i,j)=0.0d0
 334           Gcopy(i,j)=0.0d0
 335           do k=1,dimen
 336             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 337             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 338             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 339           enddo
 340         enddo
 341       enddo
 342       if (lprn) then
 343         write (iout,*) "Comparison of original and restored G"
 344         do i=1,dimen
 345           do j=1,dimen
 346             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
 347      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 348           enddo
 349         enddo
 350       endif
 351       return
 352       end
 353 c-------------------------------------------------------------------------------
 354       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
 355       implicit real*8 (a-h,o-z)
 356       include 'DIMENSIONS'
 357       include 'COMMON.IOUNITS'
 358       double precision A(LM2,LM3),B(LM2)
 359       KA=1
 360       KC=6
 361     1 KB=MIN0(KC,NC)
 362       WRITE(IOUT,600) (I,I=KA,KB)
 363       WRITE(IOUT,601) (B(I),I=KA,KB)
 364       WRITE(IOUT,602)
 365     2 N=0
 366       DO 3  I=1,NR
 367       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 368       N=N+1
 369       IF(N.LT.10) GO TO 3
 370       WRITE(IOUT,602)
 371       N=0
 372     3 CONTINUE
 373     4 IF (KB.EQ.NC) RETURN
 374       KA=KC+1
 375       KC=KC+6
 376       GO TO 1
 377   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 378   601 FORMAT (/5X,10(1PE11.4))
 379   602 FORMAT (2H  )
 380   603 FORMAT (I5,10F11.5)
 381   604 FORMAT (1H1)
 382       END
 383 c-------------------------------------------------------------------------------
 384       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
 385       implicit real*8 (a-h,o-z)
 386       include 'DIMENSIONS'
 387       include 'COMMON.IOUNITS'
 388       double precision A(LM2,LM3)
 389       KA=1
 390       KC=6
 391     1 KB=MIN0(KC,NC)
 392       WRITE(IOUT,600) (I,I=KA,KB)
 393       WRITE(IOUT,602)
 394     2 N=0
 395       DO 3  I=1,NR
 396       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 397       N=N+1
 398       IF(N.LT.10) GO TO 3
 399       WRITE(IOUT,602)
 400       N=0
 401     3 CONTINUE
 402     4 IF (KB.EQ.NC) RETURN
 403       KA=KC+1
 404       KC=KC+6
 405       GO TO 1
 406   600 FORMAT (//5x,9I11)
 407   602 FORMAT (2H  )
 408   603 FORMAT (I5,10F11.3)
 409   604 FORMAT (1H1)
 410       END
 411 c-------------------------------------------------------------------------------
 412       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
 413       implicit real*8 (a-h,o-z)
 414       include 'DIMENSIONS'
 415       include 'COMMON.IOUNITS'
 416       double precision A(LM2,LM3)
 417       KA=1
 418       KC=21
 419     1 KB=MIN0(KC,NC)
 420       WRITE(IOUT,600) (I,I=KA,KB)
 421       WRITE(IOUT,602)
 422     2 N=0
 423       DO 3  I=1,NR
 424       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 425       N=N+1
 426       IF(N.LT.3) GO TO 3
 427       WRITE(IOUT,602)
 428       N=0
 429     3 CONTINUE
 430     4 IF (KB.EQ.NC) RETURN
 431       KA=KC+1
 432       KC=KC+21
 433       GO TO 1
 434   600 FORMAT (//5x,7(3I5,2x))
 435   602 FORMAT (2H  )
 436   603 FORMAT (I5,7(3F5.1,2x))
 437   604 FORMAT (1H1)
 438       END
 439 c-------------------------------------------------------------------------------
 440       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
 441       implicit real*8 (a-h,o-z)
 442       include 'DIMENSIONS'
 443       include 'COMMON.IOUNITS'
 444       double precision A(LM2,LM3)
 445       KA=1
 446       KC=12
 447     1 KB=MIN0(KC,NC)
 448       WRITE(IOUT,600) (I,I=KA,KB)
 449       WRITE(IOUT,602)
 450     2 N=0
 451       DO 3  I=1,NR
 452       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 453       N=N+1
 454       IF(N.LT.3) GO TO 3
 455       WRITE(IOUT,602)
 456       N=0
 457     3 CONTINUE
 458     4 IF (KB.EQ.NC) RETURN
 459       KA=KC+1
 460       KC=KC+12
 461       GO TO 1
 462   600 FORMAT (//5x,4(3I9,2x))
 463   602 FORMAT (2H  )
 464   603 FORMAT (I5,4(3F9.3,2x))
 465   604 FORMAT (1H1)
 466       END
 467 c---------------------------------------------------------------------------
 468       SUBROUTINE ginv_mult(z,d_a_tmp)
 469       implicit real*8 (a-h,o-z)
 470       include 'DIMENSIONS'
 471 #ifdef MPI
 472       include 'mpif.h'
 473       integer ierr
 474 #endif
 475       include 'COMMON.SETUP'
 476       include 'COMMON.TIME1'
 477       include 'COMMON.MD'
 478       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 479      &,time01,zcopy(dimen3)
 480 #ifdef MPI
 481       if (nfgtasks.gt.1) then
 482         if (fg_rank.eq.0) then
 483 c The matching BROADCAST for fg processors is called in ERGASTULUM
 484           time00=MPI_Wtime()
 485           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 486           time_Bcast=time_Bcast+MPI_Wtime()-time00
 487 c          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 488         endif
 489 c        write (2,*) "time00",time00
 490 c        write (2,*) "Before Scatterv"
 491 c        call flush(2)
 492 c        write (2,*) "Whole z (for FG master)"
 493 c        do i=1,dimen
 494 c          write (2,*) i,z(i)
 495 c        enddo
 496 c        call MPI_Barrier(FG_COMM,IERROR)
 497         time00=MPI_Wtime()
 498         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 499      &    MPI_DOUBLE_PRECISION,
 500      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 501 c        write (2,*) "My chunk of z"
 502         do i=1,3*my_ng_count
 503           z(i)=zcopy(i)
 504 c          write (2,*) i,z(i)
 505         enddo
 506 c        write (2,*) "After SCATTERV"
 507 c        call flush(2)
 508 c        write (2,*) "MPI_Wtime",MPI_Wtime()
 509         time_scatter=time_scatter+MPI_Wtime()-time00
 510 #ifdef TIMING
 511         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 512 #endif
 513 c        write (2,*) "time_scatter",time_scatter
 514 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 515 c     &    my_ng_count
 516 c        call flush(2)
 517         time01=MPI_Wtime()
 518         do k=0,2
 519           do i=1,dimen
 520             ind=(i-1)*3+k+1
 521             temp(ind)=0.0d0
 522             do j=1,my_ng_count
 523 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 524 c     &         Ginv(i,j),z((j-1)*3+k+1),
 525 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 526 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 527               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
 528             enddo
 529           enddo
 530         enddo
 531         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 532 c        write (2,*) "Before REDUCE"
 533 c        call flush(2)
 534 c        write (2,*) "z before reduce"
 535 c        do i=1,dimen
 536 c          write (2,*) i,temp(i)
 537 c        enddo
 538         time00=MPI_Wtime()
 539         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 540      &      MPI_SUM,king,FG_COMM,IERR)
 541         time_reduce=time_reduce+MPI_Wtime()-time00
 542 c        write (2,*) "After REDUCE"
 543 c        call flush(2)
 544       else
 545 #endif
 546 #ifdef TIMING
 547         time01=MPI_Wtime()
 548 #endif
 549         do k=0,2
 550           do i=1,dimen
 551             ind=(i-1)*3+k+1
 552             d_a_tmp(ind)=0.0d0
 553             do j=1,dimen
 554 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 555 c              call flush(2)
 556 c     &         Ginv(i,j),z((j-1)*3+k+1),
 557 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 558               d_a_tmp(ind)=d_a_tmp(ind)
 559      &                         +Ginv(j,i)*z((j-1)*3+k+1)
 560 c              d_a_tmp(ind)=d_a_tmp(ind)
 561 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 562             enddo
 563           enddo
 564         enddo
 565 #ifdef TIMING
 566         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 567 #endif
 568 #ifdef MPI
 569       endif
 570 #endif
 571       return
 572       end
 573 c---------------------------------------------------------------------------
 574 #ifdef GINV_MULT
 575       SUBROUTINE ginv_mult_test(z,d_a_tmp)
 576       include 'DIMENSIONS'
 577       integer dimen
 578 c      include 'COMMON.MD'
 579       double precision z(dimen),d_a_tmp(dimen)
 580       double precision ztmp(dimen/3),dtmp(dimen/3)
 581
 582 c      do i=1,dimen
 583 c        d_a_tmp(i)=0.0d0
 584 c        do j=1,dimen
 585 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 586 c        enddo
 587 c      enddo
 588 c
 589 c      return
 590
 591 !ibm* unroll(3)
 592       do k=0,2
 593        do j=1,dimen/3
 594         ztmp(j)=z((j-1)*3+k+1)
 595        enddo
 596
 597        call alignx(16,ztmp(1))
 598        call alignx(16,dtmp(1))
 599        call alignx(16,Ginv(1,1))
 600
 601        do i=1,dimen/3
 602         dtmp(i)=0.0d0
 603         do j=1,dimen/3
 604            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 605         enddo
 606        enddo
 607        do i=1,dimen/3
 608         ind=(i-1)*3+k+1
 609         d_a_tmp(ind)=dtmp(i)
 610        enddo
 611       enddo
 612       return
 613       end
 614 #endif
 615 c---------------------------------------------------------------------------
 616       SUBROUTINE fricmat_mult(z,d_a_tmp)
 617       include 'DIMENSIONS'
 618 #ifdef MPI
 619       include 'mpif.h'
 620       integer IERROR
 621 #endif
 622       include 'COMMON.MD'
 623       include 'COMMON.IOUNITS'
 624       include 'COMMON.SETUP'
 625       include 'COMMON.TIME1'
 626 #ifndef LANG0
 627       include 'COMMON.LANGEVIN'
 628 #else
 629       include 'COMMON.LANGEVIN.lang0'
 630 #endif
 631       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 632      &,time01,zcopy(dimen3)
 633 #ifdef MPI
 634       if (nfgtasks.gt.1) then
 635         if (fg_rank.eq.0) then
 636 c The matching BROADCAST for fg processors is called in ERGASTULUM
 637           time00=MPI_Wtime()
 638           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
 639           time_Bcast=time_Bcast+MPI_Wtime()-time00
 640 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
 641         endif
 642 c        call MPI_Barrier(FG_COMM,IERROR)
 643         time00=MPI_Wtime()
 644         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 645      &    MPI_DOUBLE_PRECISION,
 646      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 647 c        write (2,*) "My chunk of z"
 648         do i=1,3*my_ng_count
 649           z(i)=zcopy(i)
 650 c          write (2,*) i,z(i)
 651         enddo
 652         time_scatter=time_scatter+MPI_Wtime()-time00
 653 #ifdef TIMING
 654         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
 655 #endif
 656         time01=MPI_Wtime()
 657         do k=0,2
 658           do i=1,dimen
 659             ind=(i-1)*3+k+1
 660             temp(ind)=0.0d0
 661             do j=1,my_ng_count
 662               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
 663             enddo
 664           enddo
 665         enddo
 666         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 667 c        write (2,*) "Before REDUCE"
 668 c        write (2,*) "d_a_tmp before reduce"
 669 c        do i=1,dimen3
 670 c          write (2,*) i,temp(i)
 671 c        enddo
 672 c        call flush(2)
 673         time00=MPI_Wtime()
 674         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 675      &      MPI_SUM,king,FG_COMM,IERR)
 676         time_reduce=time_reduce+MPI_Wtime()-time00
 677 c        write (2,*) "After REDUCE"
 678 c        call flush(2)
 679       else
 680 #endif
 681 #ifdef TIMING
 682         time01=MPI_Wtime()
 683 #endif
 684         do k=0,2
 685          do i=1,dimen
 686           ind=(i-1)*3+k+1
 687           d_a_tmp(ind)=0.0d0
 688           do j=1,dimen
 689              d_a_tmp(ind)=d_a_tmp(ind)
 690      &                           -fricmat(j,i)*z((j-1)*3+k+1)
 691           enddo
 692          enddo
 693         enddo
 694 #ifdef TIMING
 695         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 696 #endif
 697 #ifdef MPI
 698       endif
 699 #endif
 700 c      write (iout,*) "Vector d_a"
 701 c      do i=1,dimen3
 702 c        write (2,*) i,d_a_tmp(i)
 703 c      enddo
 704       return
 705       end