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