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-------------------------------------------------------------------------
14 include 'COMMON.CHAIN'
15 include 'COMMON.DERIV'
17 include 'COMMON.LOCAL'
18 include 'COMMON.INTERACT'
21 include 'COMMON.LAGRANGE.5diag'
23 include 'COMMON.LAGRANGE'
27 include 'COMMON.LANGEVIN.lang0.5diag'
29 include 'COMMON.LANGEVIN.lang0'
32 include 'COMMON.LANGEVIN'
34 include 'COMMON.IOUNITS'
35 include 'COMMON.CONTROL'
37 include 'COMMON.TIME1'
40 double precision zapas(MAXRES6),muca_factor
41 logical lprn /.false./
43 common /cipiszcze/ itime
45 double precision rs(maxres2_chain),xsolv(maxres2_chain),ip4
46 double precision aaux(3)
47 integer nind,innt,inct,inct_prev,ichain,n,mark
49 double precision rscheck(maxres2_chain),rsold(maxres2_chain)
60 write (iout,*) "Potential forces backbone"
62 write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
64 write (iout,*) "Potential forces sidechain"
66 ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
67 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
72 innt=iposd_chain(ichain)
75 do i=chain_border(1,ichain),chain_border(2,ichain)
76 if (itype(i).eq.10)then
77 rs(ind)=-gcart(j,i)-gxcart(j,i)
81 rs(ind+1)=-gxcart(j,i)
90 & "RHS of the 5-diag equations system, chain",ichain," j",j
92 write(iout,'(i5,f10.5)') i,rs(i)
95 call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv)
97 write (iout,*) "Solution of the 5-diagonal equations system"
99 write (iout,'(i5,f10.5)') i,xsolv(i)
104 call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),
107 write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),
108 & "ratio",rscheck(i)/rsold(i)
114 do i=chain_border(1,ichain),chain_border(2,ichain)
115 if (itype(i).eq.10) then
120 d_a(j,i+nres)=xsolv(ind+1)
127 write (iout,*) "Acceleration in CA and SC oordinates"
129 write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
132 write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
135 C Conevert d_a to virtual-bon-vector basis
138 c write (iout,*) "WLOS"
143 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
145 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
149 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
150 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
159 c write (iout,*) "Shifting accelerations"
167 c write (iout,*) "ichain",chain_border1(1,ichain)-1,
168 c & chain_border1(1,ichain)
169 d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
170 d_a(:,chain_border1(1,ichain))=0.0d0
172 c write (iout,*) "Adding accelerations"
174 c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
175 c & chain_border(2,ichain-1)
176 d_a(:,chain_border1(1,ichain)-1)=
177 & d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
178 d_a(:,chain_border(2,ichain-1))=0.0d0
187 innt=chain_border(1,ichain)
188 inct=chain_border(2,ichain)
190 d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
194 if (itype(i).ne.10) then
196 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
205 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
211 write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
213 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
216 write (iout,'(i3,3f10.5,3x,3f10.5)')
217 & i,(d_a(j,i+nres),j=1,3)
226 write (iout,*) "Potential forces backbone"
229 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
230 & i,(-gcart(j,i),j=1,3)
233 zapas(ind)=-gcart(j,i)
236 if (lprn) write (iout,*) "Potential forces sidechain"
238 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
239 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
240 & i,(-gcart(j,i),j=1,3)
243 zapas(ind)=-gxcart(j,i)
248 call ginv_mult(zapas,d_a_work)
257 d_a(j,i)=d_a_work(ind)
261 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
264 d_a(j,i+nres)=d_a_work(ind)
271 if(mucadyn.gt.0) call muca_update(potE)
272 factor=muca_factor(potE)*t_bath*Rb
274 cd print *,'lmuca ',factor,potE
276 d_a(j,0)=d_a(j,0)*factor
280 d_a(j,i)=d_a(j,i)*factor
285 d_a(j,i+nres)=d_a(j,i+nres)*factor
292 write(iout,*) 'acceleration 3D'
293 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
295 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
298 write (iout,'(i3,3f10.5,3x,3f10.5)')
299 & i+nres,(d_a(j,i+nres),j=1,3)
304 time_lagrangian=time_lagrangian+MPI_Wtime()-time00
308 c------------------------------------------------------------------
309 subroutine setup_MD_matrices
315 double precision time00
317 include 'COMMON.CONTROL'
318 include 'COMMON.SETUP'
320 include 'COMMON.CHAIN'
321 include 'COMMON.DERIV'
323 include 'COMMON.LOCAL'
324 include 'COMMON.INTERACT'
327 include 'COMMON.LAGRANGE.5diag'
329 include 'COMMON.LAGRANGE'
332 include 'COMMON.LANGEVIN'
335 include 'COMMON.LANGEVIN.lang0.5diag'
337 include 'COMMON.LANGEVIN.lang0'
340 include 'COMMON.IOUNITS'
341 include 'COMMON.TIME1'
342 integer i,j,k,m,m1,ind,ind1,ii,iti,ii1,jj
343 double precision coeff
344 logical lprn /.false./
346 double precision dtdi,massvec(maxres2)
348 integer ichain,innt,inct,nind,mark,n
351 double precision Gcopy(maxres2,maxres2),Ghalf(mmaxres2),
353 double precision work(8*maxres6)
354 integer iwork(maxres6)
355 common /przechowalnia/ Gcopy,Ghalf
358 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
359 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
361 c Determine the number of degrees of freedom (dimen) and the number of
367 dimen=dimen+chain_length(ichain)
368 dimen1=dimen1+2*chain_length(ichain)-1
369 dimenp=dimenp+chain_length(ichain)-1
371 write (iout,*) "Number of Calphas",dimen
372 write (iout,*) "Number of sidechains",nside
373 write (iout,*) "Number of peptide groups",dimenp
374 dimen=dimen+nside ! number of centers
375 dimen3=3*dimen ! degrees of freedom
376 write (iout,*) "Number of centers",dimen
377 write (iout,*) "Degrees of freedom:",dimen3
381 iposd_chain(ichain)=ind
382 innt=chain_border(1,ichain)
383 inct=chain_border(2,ichain)
385 if (iabs(itype(innt)).eq.10) then
386 DM(ind)=DM(ind)+msc(10)
390 DM(ind)=DM(ind)+isc(iabs(itype(innt)))
391 DM(ind+1)=msc(iabs(itype(innt)))+isc(iabs(itype(innt)))
395 write (iout,*) "ind",ind," nind",nind
397 ! if (iabs(itype(i)).eq.ntyp1) cycle
399 if (iabs(itype(i)).eq.10) then
400 if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
404 DM(ind)=DM(ind)+isc(iabs(itype(i)))
405 DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
409 write (iout,*) "i",i," ind",ind," nind",nind
411 if (inct.gt.innt) then
413 if (iabs(itype(inct)).eq.10) then
414 DM(ind)=DM(ind)+msc(10)
418 DM(ind)=DM(ind)+isc(iabs(itype(inct)))
419 DM(ind+1)=msc(iabs(itype(inct)))+isc(iabs(itype(inct)))
424 write (iout,*) "ind",ind," nind",nind
425 dimen_chain(ichain)=nind
429 ind=iposd_chain(ichain)
430 innt=chain_border(1,ichain)
431 inct=chain_border(2,ichain)
433 if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
434 DU1(ind)=-isc(iabs(itype(i)))
445 ind=iposd_chain(ichain)
446 innt=chain_border(1,ichain)
447 inct=chain_border(2,ichain)
449 ! if (iabs(itype(i)).eq.ntyp1) cycle
450 c write (iout,*) "i",i," itype",itype(i),ntyp1
451 if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
466 write (iout,*)"The upper part of the five-diagonal inertia matrix"
469 if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
470 n=dimen_chain(ichain)
471 innt=iposd_chain(ichain)
472 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
475 if (i.lt.inct-1) then
476 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
477 else if (i.eq.inct-1) then
478 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
480 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
484 call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),
485 & DU2(innt:inct-1), MARK)
489 & "ERROR: the inertia matrix is not positive definite for chain",
492 call MPI_Finalize(ierr)
495 else if (mark.eq.0) then
497 & "ERROR: the inertia matrix is singular for chain",ichain
499 call MPI_Finalize(ierr)
501 else if (mark.eq.1) then
503 write (iout,*) "The transformed five-diagonal inertia matrix"
504 write (iout,'(a,i5)') 'Chain',ichain
506 if (i.lt.inct-1) then
507 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
508 else if (i.eq.inct-1) then
509 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
511 write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
517 ! Diagonalization of the pentadiagonal matrix
522 dimen=(nct-nnt+1)+nside
523 dimen1=(nct-nnt)+(nct-nnt+1)
525 write (iout,*) "Degrees_of_freedom",dimen3
527 if (nfgtasks.gt.1) then
529 call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
530 time_Bcast=time_Bcast+MPI_Wtime()-time00
531 call int_bounds(dimen,igmult_start,igmult_end)
532 igmult_start=igmult_start-1
533 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
534 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
535 my_ng_count=igmult_end-igmult_start
536 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
537 & MPI_INTEGER,FG_COMM,IERROR)
538 write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
539 & ' absolute rank',myrank,' igmult_start',igmult_start,
540 & ' igmult_end',igmult_end,' count',my_ng_count
541 write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
542 write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
552 c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
553 c Zeroing out A and fricmat
559 c Diagonal elements of the dC part of A and the respective friction coefficients
571 c Off-diagonal elements of the dC part of A
578 c Diagonal elements of the dX part of A and the respective friction coefficients
588 massvec(ii)=msc(iabs(iti))
589 if (iti.ne.10 .and. iti.ne.ntyp1) then
593 Gmat(ii1,ii1)=ISC(iabs(iti))
596 c Off-diagonal elements of the dX part of A
610 write (iout,*) "Vector massvec"
612 write (iout,*) i,massvec(i)
614 write (iout,'(//a)') "A"
615 call matout(dimen,dimen1,maxres2,maxres2,A)
618 c Calculate the G matrix (store in Gmat)
623 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
625 Gmat(k,i)=Gmat(k,i)+dtdi
630 write (iout,'(//a)') "Gmat"
631 call matout(dimen,dimen,maxres2,maxres2,Gmat)
640 c Invert the G matrix
641 call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
643 write (iout,'(//a)') "Ginv"
644 call matout(dimen,dimen,maxres2,maxres2,Ginv)
647 if (nfgtasks.gt.1) then
648 myginv_ng_count=maxres2*my_ng_count
649 call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
650 & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
651 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
652 & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
653 write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
654 write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
656 c call MPI_Scatterv(ginv(1,1),nginv_counts(0),
657 c & nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
658 c & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
659 c call MPI_Barrier(FG_COMM,IERR)
661 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
662 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
663 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
665 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
672 c write (iout,*) "Master's chunk of ginv"
673 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
677 write (iout,*) "The G matrix is singular."
680 c Compute G**(-1/2) and G**(1/2)
688 call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
692 & "Eigenvectors and eigenvalues of the G matrix"
693 call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
696 sqreig(i)=dsqrt(Geigen(i))
704 Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
705 Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
706 Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
711 write (iout,*) "Comparison of original and restored G"
714 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
715 & Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
722 c-------------------------------------------------------------------------------
723 SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
726 include 'COMMON.IOUNITS'
727 double precision A(LM2,LM3),B(LM2)
728 integer nc,nr,lm2,lm3,ka,kb,kc,n,i,j
732 WRITE(IOUT,600) (I,I=KA,KB)
733 WRITE(IOUT,601) (B(I),I=KA,KB)
737 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
743 4 IF (KB.EQ.NC) RETURN
747 600 FORMAT (// 9H ROOT NO.,I4,9I11)
748 601 FORMAT (/5X,10(1PE11.4))
750 603 FORMAT (I5,10F11.5)
753 c-------------------------------------------------------------------------------
754 SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
757 include 'COMMON.IOUNITS'
758 double precision A(LM2,LM3)
759 integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
763 WRITE(IOUT,600) (I,I=KA,KB)
767 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
773 4 IF (KB.EQ.NC) RETURN
777 600 FORMAT (//5x,9I11)
779 603 FORMAT (I5,10F11.3)
782 c-------------------------------------------------------------------------------
783 SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
786 include 'COMMON.IOUNITS'
787 double precision A(LM2,LM3)
788 integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
792 WRITE(IOUT,600) (I,I=KA,KB)
796 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
802 4 IF (KB.EQ.NC) RETURN
806 600 FORMAT (//5x,7(3I5,2x))
808 603 FORMAT (I5,7(3F5.1,2x))
811 c-------------------------------------------------------------------------------
812 SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
815 include 'COMMON.IOUNITS'
816 double precision A(LM2,LM3)
817 integer nc,nr,lm2,lm3,ka,kb,kc,i,j,n
821 WRITE(IOUT,600) (I,I=KA,KB)
825 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
831 4 IF (KB.EQ.NC) RETURN
835 600 FORMAT (//5x,4(3I9,2x))
837 603 FORMAT (I5,4(3F9.3,2x))
840 c-----------------------------------------------------------------------------
841 SUBROUTINE MATOUTR(N,A)
842 c Prints the lower fragment of a symmetric matix
845 double precision a(n*(n+1)/2)
846 integer i,j,k,nlim,jlim,jlim1
847 CHARACTER*6 LINE6 / '------' /
848 CHARACTER*12 LINE12 / '------------' /
849 double precision B(10)
850 include 'COMMON.IOUNITS'
853 WRITE (IOUT,1000) (K,K=I,NLIM)
854 WRITE (IOUT,1020) LINE6,(LINE12,K=I,NLIM)
855 1000 FORMAT (/7X,10(I5,2X))
856 1020 FORMAT (A6,10A7)
861 3 B(K-I+1)=A(J*(J-1)/2+K)
862 WRITE (IOUT,1010) J,(B(K),K=1,JLIM1)
865 1010 FORMAT (I3,3X,10(F7.2))
869 c---------------------------------------------------------------------------
870 subroutine fivediagmult(n,DM,DU1,DU2,x,y)
873 double precision DM(n),DU1(n),DU2(n),x(n),y(n)
875 y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
876 y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
878 y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i)
879 & +DU1(i)*x(i+1)+DU2(i)*x(i+2)
881 y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1)
883 y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
886 c---------------------------------------------------------------------------
887 subroutine fivediaginv_mult(ndim,forces,d_a_vec)
890 include 'COMMON.CHAIN'
891 include 'COMMON.IOUNITS'
892 include 'COMMON.LAGRANGE.5diag'
893 include 'COMMON.INTERACT'
896 double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
897 & xsolv(ndim),d_a_vec(6*nres)
898 integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
901 Compute accelerations in Calpha and SC
903 n=dimen_chain(ichain)
904 iposc=iposd_chain(ichain)
905 innt=chain_border(1,ichain)
906 inct=chain_border(2,ichain)
908 rs(i-iposc+1)=forces(3*(i-1)+j)
911 write (iout,*) "j",j," chain",ichain
913 write (iout,'(f10.5)') (rs(i),i=1,n)
915 call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
917 write (iout,*) "xsolv"
918 write (iout,'(f10.5)') (xsolv(i),i=1,n)
922 if (itype(i).eq.10)then
923 accel(j,i)=xsolv(ind)
926 accel(j,i)=xsolv(ind)
927 accel(j,i+nres)=xsolv(ind+1)
933 C Convert d_a to virtual-bon-vector basis
935 write (iout,*) "accel in CA-SC basis"
937 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
938 & (accel(j,i+nres),j=1,3)
940 write (iout,*) "nnt",nnt
943 accel(:,0)=accel(:,1)
946 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
948 accel(j,i)=accel(j,i+1)-accel(j,i)
952 accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
953 accel(j,i)=accel(j,i+1)-accel(j,i)
959 accel(:,2*nres)=0.0d0
961 accel(:,0)=accel(:,1)
965 accel(:,chain_border1(1,ichain)-1)=
966 & accel(:,chain_border1(1,ichain))
967 accel(:,chain_border1(1,ichain))=0.0d0
970 accel(:,chain_border1(1,ichain)-1)=
971 & accel(:,chain_border1(1,ichain)-1)
972 & +accel(:,chain_border(2,ichain-1))
973 accel(:,chain_border(2,ichain-1))=0.0d0
976 write (iout,*) "accel in fivediaginv_mult: 1"
978 write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
982 d_a_vec(j)=accel(j,0)
987 d_a_vec(ind+j)=accel(j,i)
992 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
994 d_a_vec(ind+j)=accel(j,i+nres)
1000 write (iout,*) "d_a_vec"
1001 write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
1006 c---------------------------------------------------------------------------
1007 SUBROUTINE ginv_mult(z,d_a_tmp)
1009 include 'DIMENSIONS'
1014 include 'COMMON.SETUP'
1015 include 'COMMON.TIME1'
1017 include 'COMMON.LAGRANGE'
1018 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1019 &,time01,zcopy(dimen3)
1022 if (nfgtasks.gt.1) then
1023 if (fg_rank.eq.0) then
1024 c The matching BROADCAST for fg processors is called in ERGASTULUM
1026 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
1027 time_Bcast=time_Bcast+MPI_Wtime()-time00
1028 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
1030 c write (2,*) "time00",time00
1031 c write (2,*) "Before Scatterv"
1033 c write (2,*) "Whole z (for FG master)"
1035 c write (2,*) i,z(i)
1037 c call MPI_Barrier(FG_COMM,IERROR)
1039 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1040 & MPI_DOUBLE_PRECISION,
1041 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1042 c write (2,*) "My chunk of z"
1043 do i=1,3*my_ng_count
1045 c write (2,*) i,z(i)
1047 c write (2,*) "After SCATTERV"
1049 c write (2,*) "MPI_Wtime",MPI_Wtime()
1050 time_scatter=time_scatter+MPI_Wtime()-time00
1052 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1054 c write (2,*) "time_scatter",time_scatter
1055 c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1064 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1065 c & Ginv(i,j),z((j-1)*3+k+1),
1066 c & Ginv(i,j)*z((j-1)*3+k+1)
1067 c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1068 temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
1072 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1073 c write (2,*) "Before REDUCE"
1075 c write (2,*) "z before reduce"
1077 c write (2,*) i,temp(i)
1080 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1081 & MPI_SUM,king,FG_COMM,IERR)
1082 time_reduce=time_reduce+MPI_Wtime()-time00
1083 c write (2,*) "After REDUCE"
1095 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1097 c & Ginv(i,j),z((j-1)*3+k+1),
1098 c & Ginv(i,j)*z((j-1)*3+k+1)
1099 d_a_tmp(ind)=d_a_tmp(ind)
1100 & +Ginv(j,i)*z((j-1)*3+k+1)
1101 c d_a_tmp(ind)=d_a_tmp(ind)
1102 c & +Ginv(i,j)*z((j-1)*3+k+1)
1107 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1114 c---------------------------------------------------------------------------
1116 SUBROUTINE ginv_mult_test(z,d_a_tmp)
1118 include 'DIMENSIONS'
1119 include 'COMMON.LAGRANGE'
1120 double precision z(dimen),d_a_tmp(dimen)
1121 double precision ztmp(dimen/3),dtmp(dimen/3)
1126 c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1135 ztmp(j)=z((j-1)*3+k+1)
1138 call alignx(16,ztmp(1))
1139 call alignx(16,dtmp(1))
1140 call alignx(16,Ginv(1,1))
1145 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1150 d_a_tmp(ind)=dtmp(i)
1156 c---------------------------------------------------------------------------
1157 SUBROUTINE fricmat_mult(z,d_a_tmp)
1158 include 'DIMENSIONS'
1164 include 'COMMON.LAGRANGE'
1165 include 'COMMON.IOUNITS'
1166 include 'COMMON.SETUP'
1167 include 'COMMON.TIME1'
1169 include 'COMMON.LANGEVIN'
1171 include 'COMMON.LANGEVIN.lang0'
1173 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1174 &,time01,zcopy(dimen3)
1176 if (nfgtasks.gt.1) then
1177 if (fg_rank.eq.0) then
1178 c The matching BROADCAST for fg processors is called in ERGASTULUM
1180 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1181 time_Bcast=time_Bcast+MPI_Wtime()-time00
1182 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1184 c call MPI_Barrier(FG_COMM,IERROR)
1186 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1187 & MPI_DOUBLE_PRECISION,
1188 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1189 c write (2,*) "My chunk of z"
1190 do i=1,3*my_ng_count
1192 c write (2,*) i,z(i)
1194 time_scatter=time_scatter+MPI_Wtime()-time00
1196 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1204 temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
1208 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1209 c write (2,*) "Before REDUCE"
1210 c write (2,*) "d_a_tmp before reduce"
1212 c write (2,*) i,temp(i)
1216 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1217 & MPI_SUM,king,FG_COMM,IERR)
1218 time_reduce=time_reduce+MPI_Wtime()-time00
1219 c write (2,*) "After REDUCE"
1231 d_a_tmp(ind)=d_a_tmp(ind)
1232 & -fricmat(j,i)*z((j-1)*3+k+1)
1237 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1242 c write (iout,*) "Vector d_a"
1244 c write (2,*) i,d_a_tmp(i)