Merge branch 'homology' of mmka.chem.univ.gda.pl:unres into homology
[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.21) 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.21) 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(21)=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
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      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
501 c        write (2,*) "My chunk of z"
502 c        do i=1,3*my_ng_count
503 c          write (2,*) i,z(i)
504 c        enddo
505 c        write (2,*) "After SCATTERV"
506 c        call flush(2)
507 c        write (2,*) "MPI_Wtime",MPI_Wtime()
508         time_scatter=time_scatter+MPI_Wtime()-time00
509 #ifdef TIMING
510         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
511 #endif
512 c        write (2,*) "time_scatter",time_scatter
513 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
514 c     &    my_ng_count
515 c        call flush(2)
516         time01=MPI_Wtime()
517         do k=0,2
518           do i=1,dimen
519             ind=(i-1)*3+k+1
520             temp(ind)=0.0d0
521             do j=1,my_ng_count
522 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
523 c     &         Ginv(i,j),z((j-1)*3+k+1),
524 c     &          Ginv(i,j)*z((j-1)*3+k+1)
525 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
526               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
527             enddo
528           enddo 
529         enddo
530         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
531 c        write (2,*) "Before REDUCE"
532 c        call flush(2)
533 c        write (2,*) "z before reduce"
534 c        do i=1,dimen
535 c          write (2,*) i,temp(i)
536 c        enddo
537         time00=MPI_Wtime()
538         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
539      &      MPI_SUM,king,FG_COMM,IERR)
540         time_reduce=time_reduce+MPI_Wtime()-time00
541 c        write (2,*) "After REDUCE"
542 c        call flush(2)
543       else
544 #endif
545 #ifdef TIMING
546         time01=MPI_Wtime()
547 #endif
548         do k=0,2
549           do i=1,dimen
550             ind=(i-1)*3+k+1
551             d_a_tmp(ind)=0.0d0
552             do j=1,dimen
553 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
554 c              call flush(2)
555 c     &         Ginv(i,j),z((j-1)*3+k+1),
556 c     &          Ginv(i,j)*z((j-1)*3+k+1)
557               d_a_tmp(ind)=d_a_tmp(ind)
558      &                         +Ginv(j,i)*z((j-1)*3+k+1)
559 c              d_a_tmp(ind)=d_a_tmp(ind)
560 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
561             enddo
562           enddo 
563         enddo
564 #ifdef TIMING
565         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
566 #endif
567 #ifdef MPI
568       endif
569 #endif
570       return
571       end
572 c---------------------------------------------------------------------------
573 #ifdef GINV_MULT
574       SUBROUTINE ginv_mult_test(z,d_a_tmp)
575       include 'DIMENSIONS'
576       integer dimen
577 c      include 'COMMON.MD'
578       double precision z(dimen),d_a_tmp(dimen)
579       double precision ztmp(dimen/3),dtmp(dimen/3)
580
581 c      do i=1,dimen
582 c        d_a_tmp(i)=0.0d0
583 c        do j=1,dimen
584 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
585 c        enddo
586 c      enddo
587 c
588 c      return
589
590 !ibm* unroll(3)
591       do k=0,2
592        do j=1,dimen/3
593         ztmp(j)=z((j-1)*3+k+1)
594        enddo
595
596        call alignx(16,ztmp(1))
597        call alignx(16,dtmp(1))
598        call alignx(16,Ginv(1,1)) 
599
600        do i=1,dimen/3
601         dtmp(i)=0.0d0
602         do j=1,dimen/3
603            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
604         enddo
605        enddo
606        do i=1,dimen/3
607         ind=(i-1)*3+k+1
608         d_a_tmp(ind)=dtmp(i)
609        enddo 
610       enddo
611       return
612       end
613 #endif
614 c---------------------------------------------------------------------------
615       SUBROUTINE fricmat_mult(z,d_a_tmp)
616       include 'DIMENSIONS'
617 #ifdef MPI
618       include 'mpif.h'
619       integer IERROR
620 #endif
621       include 'COMMON.MD'
622       include 'COMMON.IOUNITS'
623       include 'COMMON.SETUP'
624       include 'COMMON.TIME1'
625 #ifndef LANG0
626       include 'COMMON.LANGEVIN'
627 #else
628       include 'COMMON.LANGEVIN.lang0'
629 #endif
630       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
631      &time01
632 #ifdef MPI
633       if (nfgtasks.gt.1) then
634         if (fg_rank.eq.0) then
635 c The matching BROADCAST for fg processors is called in ERGASTULUM
636           time00=MPI_Wtime()
637           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
638           time_Bcast=time_Bcast+MPI_Wtime()-time00
639 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
640         endif
641 c        call MPI_Barrier(FG_COMM,IERROR)
642         time00=MPI_Wtime()
643         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
644      &    MPI_DOUBLE_PRECISION,
645      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
646 c        write (2,*) "My chunk of z"
647 c        do i=1,3*my_ng_count
648 c          write (2,*) i,z(i)
649 c        enddo
650         time_scatter=time_scatter+MPI_Wtime()-time00
651 #ifdef TIMING
652         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
653 #endif
654         time01=MPI_Wtime()
655         do k=0,2
656           do i=1,dimen
657             ind=(i-1)*3+k+1
658             temp(ind)=0.0d0
659             do j=1,my_ng_count
660               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
661             enddo
662           enddo 
663         enddo
664         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
665 c        write (2,*) "Before REDUCE"
666 c        write (2,*) "d_a_tmp before reduce"
667 c        do i=1,dimen3
668 c          write (2,*) i,temp(i)
669 c        enddo
670 c        call flush(2)
671         time00=MPI_Wtime()
672         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
673      &      MPI_SUM,king,FG_COMM,IERR)
674         time_reduce=time_reduce+MPI_Wtime()-time00
675 c        write (2,*) "After REDUCE"
676 c        call flush(2)
677       else
678 #endif
679 #ifdef TIMING
680         time01=MPI_Wtime()
681 #endif
682         do k=0,2
683          do i=1,dimen
684           ind=(i-1)*3+k+1
685           d_a_tmp(ind)=0.0d0
686           do j=1,dimen
687              d_a_tmp(ind)=d_a_tmp(ind)
688      &                           -fricmat(j,i)*z((j-1)*3+k+1)
689           enddo
690          enddo 
691         enddo
692 #ifdef TIMING
693         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
694 #endif
695 #ifdef MPI
696       endif
697 #endif
698 c      write (iout,*) "Vector d_a"
699 c      do i=1,dimen3
700 c        write (2,*) i,d_a_tmp(i)
701 c      enddo
702       return
703       end