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