changes in wham
[unres.git] / 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