compinfo.c rewrite
[unres.git] / 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 c  Zeroing out A and fricmat
188       do i=1,dimen
189         do j=1,dimen
190           A(i,j)=0.0D0     
191         enddo    
192       enddo
193 c  Diagonal elements of the dC part of A and the respective friction coefficients
194       ind=1
195       ind1=0
196       do i=nnt,nct-1
197         ind=ind+1
198         ind1=ind1+1
199         coeff=0.25d0*IP
200         massvec(ind1)=mp
201         Gmat(ind,ind)=coeff
202         A(ind1,ind)=0.5d0
203       enddo
204       
205 c  Off-diagonal elements of the dC part of A 
206       k=3
207       do i=1,nct-nnt
208         do j=1,i
209           A(i,j)=1.0d0
210         enddo
211       enddo
212 c  Diagonal elements of the dX part of A and the respective friction coefficients
213       m=nct-nnt
214       m1=nct-nnt+1
215       ind=0
216       ind1=0
217       do i=nnt,nct
218         ind=ind+1
219         ii = ind+m
220         iti=itype(i)
221         massvec(ii)=msc(iti)
222         if (iti.ne.10) then
223           ind1=ind1+1
224           ii1= ind1+m1
225           A(ii,ii1)=1.0d0
226           Gmat(ii1,ii1)=ISC(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
331       do i=1,dimen
332         sqreig(i)=dsqrt(Geigen(i))
333         invsqreig(i)=1.d0/sqreig(i)
334       enddo
335       do i=1,dimen
336          do j=1,dimen
337             Gvectmp(i,j)=Gvec(j,i)
338          enddo
339       enddo
340
341       do i=1,dimen
342         do j=1,dimen
343           Gsqrptmp=0.0d0
344           Gsqrmtmp=0.0d0
345           Gcopytmp=0.0d0
346           do k=1,dimen
347 c             Gvec2tmp=Gvec(i,k)*Gvec(j,k)
348              Gvec2tmp=Gvec(k,i)*Gvec(k,j)
349              Gsqrptmp=Gsqrptmp+Gvec2tmp*sqreig(k)
350              Gsqrmtmp=Gsqrmtmp+Gvec2tmp*invsqreig(k)
351              Gcopytmp=Gcopytmp+Gvec2tmp*Geigen(k)
352           enddo
353           Gsqrp(i,j)=Gsqrptmp
354           Gsqrm(i,j)=Gsqrmtmp
355           Gcopy(i,j)=Gcopytmp
356         enddo
357       enddo
358
359       do i=1,dimen
360          do j=1,dimen
361             Gvec(i,j)=Gvectmp(j,i)
362          enddo
363       enddo
364
365       if (lprn) then
366         write (iout,*) "Comparison of original and restored G"
367         do i=1,dimen
368           do j=1,dimen
369             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
370      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
371           enddo
372         enddo
373       endif
374       return
375       end 
376 c-------------------------------------------------------------------------------
377       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
378       implicit real*8 (a-h,o-z)
379       include 'DIMENSIONS'
380       include 'COMMON.IOUNITS'
381       double precision A(LM2,LM3),B(LM2)
382       KA=1
383       KC=6
384     1 KB=MIN0(KC,NC)
385       WRITE(IOUT,600) (I,I=KA,KB)
386       WRITE(IOUT,601) (B(I),I=KA,KB)
387       WRITE(IOUT,602)
388     2 N=0
389       DO 3  I=1,NR
390       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
391       N=N+1
392       IF(N.LT.10) GO TO 3
393       WRITE(IOUT,602)
394       N=0
395     3 CONTINUE
396     4 IF (KB.EQ.NC) RETURN
397       KA=KC+1
398       KC=KC+6
399       GO TO 1
400   600 FORMAT (// 9H ROOT NO.,I4,9I11)
401   601 FORMAT (/5X,10(1PE11.4))
402   602 FORMAT (2H  )
403   603 FORMAT (I5,10F11.5)
404   604 FORMAT (1H1)
405       END
406 c-------------------------------------------------------------------------------
407       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
408       implicit real*8 (a-h,o-z)
409       include 'DIMENSIONS'
410       include 'COMMON.IOUNITS'
411       double precision A(LM2,LM3)
412       KA=1
413       KC=6
414     1 KB=MIN0(KC,NC)
415       WRITE(IOUT,600) (I,I=KA,KB)
416       WRITE(IOUT,602)
417     2 N=0
418       DO 3  I=1,NR
419       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
420       N=N+1
421       IF(N.LT.10) GO TO 3
422       WRITE(IOUT,602)
423       N=0
424     3 CONTINUE
425     4 IF (KB.EQ.NC) RETURN
426       KA=KC+1
427       KC=KC+6
428       GO TO 1
429   600 FORMAT (//5x,9I11)
430   602 FORMAT (2H  )
431   603 FORMAT (I5,10F11.3)
432   604 FORMAT (1H1)
433       END
434 c-------------------------------------------------------------------------------
435       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
436       implicit real*8 (a-h,o-z)
437       include 'DIMENSIONS'
438       include 'COMMON.IOUNITS'
439       double precision A(LM2,LM3)
440       KA=1
441       KC=21
442     1 KB=MIN0(KC,NC)
443       WRITE(IOUT,600) (I,I=KA,KB)
444       WRITE(IOUT,602)
445     2 N=0
446       DO 3  I=1,NR
447       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
448       N=N+1
449       IF(N.LT.3) GO TO 3
450       WRITE(IOUT,602)
451       N=0
452     3 CONTINUE
453     4 IF (KB.EQ.NC) RETURN
454       KA=KC+1
455       KC=KC+21
456       GO TO 1
457   600 FORMAT (//5x,7(3I5,2x))
458   602 FORMAT (2H  )
459   603 FORMAT (I5,7(3F5.1,2x))
460   604 FORMAT (1H1)
461       END
462 c-------------------------------------------------------------------------------
463       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
464       implicit real*8 (a-h,o-z)
465       include 'DIMENSIONS'
466       include 'COMMON.IOUNITS'
467       double precision A(LM2,LM3)
468       KA=1
469       KC=12
470     1 KB=MIN0(KC,NC)
471       WRITE(IOUT,600) (I,I=KA,KB)
472       WRITE(IOUT,602)
473     2 N=0
474       DO 3  I=1,NR
475       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
476       N=N+1
477       IF(N.LT.3) GO TO 3
478       WRITE(IOUT,602)
479       N=0
480     3 CONTINUE
481     4 IF (KB.EQ.NC) RETURN
482       KA=KC+1
483       KC=KC+12
484       GO TO 1
485   600 FORMAT (//5x,4(3I9,2x))
486   602 FORMAT (2H  )
487   603 FORMAT (I5,4(3F9.3,2x))
488   604 FORMAT (1H1)
489       END
490 c---------------------------------------------------------------------------
491       SUBROUTINE ginv_mult(z,d_a_tmp)
492       implicit real*8 (a-h,o-z)
493       include 'DIMENSIONS'
494 #ifdef MPI
495       include 'mpif.h'
496       integer ierr
497 #endif
498       include 'COMMON.SETUP'
499       include 'COMMON.TIME1'
500       include 'COMMON.MD'
501       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
502      &time01
503 #ifdef MPI
504       if (nfgtasks.gt.1) then
505         if (fg_rank.eq.0) then
506 c The matching BROADCAST for fg processors is called in ERGASTULUM
507           time00=MPI_Wtime()
508           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
509           time_Bcast=time_Bcast+MPI_Wtime()-time00
510 c          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
511         endif
512 c        write (2,*) "time00",time00
513 c        write (2,*) "Before Scatterv"
514 c        call flush(2)
515 c        write (2,*) "Whole z (for FG master)"
516 c        do i=1,dimen
517 c          write (2,*) i,z(i)
518 c        enddo
519 c        call MPI_Barrier(FG_COMM,IERROR)
520         time00=MPI_Wtime()
521         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
522      &    MPI_DOUBLE_PRECISION,
523      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
524 c        write (2,*) "My chunk of z"
525 c        do i=1,3*my_ng_count
526 c          write (2,*) i,z(i)
527 c        enddo
528 c        write (2,*) "After SCATTERV"
529 c        call flush(2)
530 c        write (2,*) "MPI_Wtime",MPI_Wtime()
531         time_scatter=time_scatter+MPI_Wtime()-time00
532 #ifdef TIMING
533         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
534 #endif
535 c        write (2,*) "time_scatter",time_scatter
536 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
537 c     &    my_ng_count
538 c        call flush(2)
539         time01=MPI_Wtime()
540         do k=0,2
541           do i=1,dimen
542             ind=(i-1)*3+k+1
543             temp(ind)=0.0d0
544             do j=1,my_ng_count
545 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
546 c     &         Ginv(i,j),z((j-1)*3+k+1),
547 c     &          Ginv(i,j)*z((j-1)*3+k+1)
548 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
549               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
550             enddo
551           enddo 
552         enddo
553         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
554 c        write (2,*) "Before REDUCE"
555 c        call flush(2)
556 c        write (2,*) "z before reduce"
557 c        do i=1,dimen
558 c          write (2,*) i,temp(i)
559 c        enddo
560         time00=MPI_Wtime()
561         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
562      &      MPI_SUM,king,FG_COMM,IERR)
563         time_reduce=time_reduce+MPI_Wtime()-time00
564 c        write (2,*) "After REDUCE"
565 c        call flush(2)
566       else
567 #endif
568 #ifdef TIMING
569         time01=MPI_Wtime()
570 #endif
571         do k=0,2
572           do i=1,dimen
573             ind=(i-1)*3+k+1
574             d_a_tmp(ind)=0.0d0
575             do j=1,dimen
576 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
577 c              call flush(2)
578 c     &         Ginv(i,j),z((j-1)*3+k+1),
579 c     &          Ginv(i,j)*z((j-1)*3+k+1)
580               d_a_tmp(ind)=d_a_tmp(ind)
581      &                         +Ginv(j,i)*z((j-1)*3+k+1)
582 c              d_a_tmp(ind)=d_a_tmp(ind)
583 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
584             enddo
585           enddo 
586         enddo
587 #ifdef TIMING
588         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
589 #endif
590 #ifdef MPI
591       endif
592 #endif
593       return
594       end
595 c---------------------------------------------------------------------------
596 #ifdef GINV_MULT
597       SUBROUTINE ginv_mult_test(z,d_a_tmp)
598       include 'DIMENSIONS'
599       integer dimen
600 c      include 'COMMON.MD'
601       double precision z(dimen),d_a_tmp(dimen)
602       double precision ztmp(dimen/3),dtmp(dimen/3)
603
604 c      do i=1,dimen
605 c        d_a_tmp(i)=0.0d0
606 c        do j=1,dimen
607 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
608 c        enddo
609 c      enddo
610 c
611 c      return
612
613 !ibm* unroll(3)
614       do k=0,2
615        do j=1,dimen/3
616         ztmp(j)=z((j-1)*3+k+1)
617        enddo
618
619        call alignx(16,ztmp(1))
620        call alignx(16,dtmp(1))
621        call alignx(16,Ginv(1,1)) 
622
623        do i=1,dimen/3
624         dtmp(i)=0.0d0
625         do j=1,dimen/3
626            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
627         enddo
628        enddo
629        do i=1,dimen/3
630         ind=(i-1)*3+k+1
631         d_a_tmp(ind)=dtmp(i)
632        enddo 
633       enddo
634       return
635       end
636 #endif
637 c---------------------------------------------------------------------------
638       SUBROUTINE fricmat_mult(z,d_a_tmp)
639       include 'DIMENSIONS'
640 #ifdef MPI
641       include 'mpif.h'
642       integer IERROR
643 #endif
644       include 'COMMON.MD'
645       include 'COMMON.IOUNITS'
646       include 'COMMON.SETUP'
647       include 'COMMON.TIME1'
648 #ifndef LANG0
649       include 'COMMON.LANGEVIN'
650 #else
651       include 'COMMON.LANGEVIN.lang0'
652 #endif
653       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
654      &time01
655 #ifdef MPI
656       if (nfgtasks.gt.1) then
657         if (fg_rank.eq.0) then
658 c The matching BROADCAST for fg processors is called in ERGASTULUM
659           time00=MPI_Wtime()
660           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
661           time_Bcast=time_Bcast+MPI_Wtime()-time00
662 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
663         endif
664 c        call MPI_Barrier(FG_COMM,IERROR)
665         time00=MPI_Wtime()
666         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
667      &    MPI_DOUBLE_PRECISION,
668      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
669 c        write (2,*) "My chunk of z"
670 c        do i=1,3*my_ng_count
671 c          write (2,*) i,z(i)
672 c        enddo
673         time_scatter=time_scatter+MPI_Wtime()-time00
674 #ifdef TIMING
675         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
676 #endif
677         time01=MPI_Wtime()
678         do k=0,2
679           do i=1,dimen
680             ind=(i-1)*3+k+1
681             temp(ind)=0.0d0
682             do j=1,my_ng_count
683               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
684             enddo
685           enddo 
686         enddo
687         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
688 c        write (2,*) "Before REDUCE"
689 c        write (2,*) "d_a_tmp before reduce"
690 c        do i=1,dimen3
691 c          write (2,*) i,temp(i)
692 c        enddo
693 c        call flush(2)
694         time00=MPI_Wtime()
695         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
696      &      MPI_SUM,king,FG_COMM,IERR)
697         time_reduce=time_reduce+MPI_Wtime()-time00
698 c        write (2,*) "After REDUCE"
699 c        call flush(2)
700       else
701 #endif
702 #ifdef TIMING
703         time01=MPI_Wtime()
704 #endif
705         do k=0,2
706          do i=1,dimen
707           ind=(i-1)*3+k+1
708           d_a_tmp(ind)=0.0d0
709           do j=1,dimen
710              d_a_tmp(ind)=d_a_tmp(ind)
711      &                           -fricmat(j,i)*z((j-1)*3+k+1)
712           enddo
713          enddo 
714         enddo
715 #ifdef TIMING
716         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
717 #endif
718 #ifdef MPI
719       endif
720 #endif
721 c      write (iout,*) "Vector d_a"
722 c      do i=1,dimen3
723 c        write (2,*) i,d_a_tmp(i)
724 c      enddo
725       return
726       end