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