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