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