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