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'
895 double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
896 & xsolv(ndim),d_a_vec(6*nres)
897 integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
899 Compute accelerations in Calpha and SC
901 n=dimen_chain(ichain)
902 iposc=iposd_chain(ichain)
903 innt=chain_border(1,ichain)
904 inct=chain_border(2,ichain)
906 rs(i)=forces(3*(i-1)+j)
908 call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
911 if (itype(i).eq.10)then
912 accel(j,i)=xsolv(ind)
915 accel(j,i)=xsolv(ind)
916 accel(j,i+nres)=xsolv(ind+1)
922 C Conevert d_a to virtual-bon-vector basis
924 write (iout,*) "accel in CA-SC basis"
926 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
927 & (accel(j,i+nres),j=1,3)
929 write (iout,*) "nnt",nnt
932 accel(:,0)=accel(:,1)
935 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
937 accel(j,i)=accel(j,i+1)-accel(j,i)
941 accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
942 accel(j,i)=accel(j,i+1)-accel(j,i)
947 accel(:,2*nres)=0.0d0
949 accel(:,0)=accel(:,1)
953 accel(:,chain_border1(1,ichain)-1)=
954 & accel(:,chain_border1(1,ichain))
955 accel(:,chain_border1(1,ichain))=0.0d0
958 accel(:,chain_border1(1,ichain)-1)=
959 & accel(:,chain_border1(1,ichain)-1)
960 & +accel(:,chain_border(2,ichain-1))
961 accel(:,chain_border(2,ichain-1))=0.0d0
964 write (iout,*) "accel in fivediaginv_mult: 1"
966 write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
970 d_a_vec(j)=accel(j,0)
975 d_a_vec(ind+j)=accel(j,i)
980 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
982 d_a_vec(ind+j)=accel(j,i+nres)
988 write (iout,*) "d_a_vec"
989 write (iout,'(3f10.5)') (d_a_vec(j),j=1,dimen3)
994 c---------------------------------------------------------------------------
995 SUBROUTINE ginv_mult(z,d_a_tmp)
1002 include 'COMMON.SETUP'
1003 include 'COMMON.TIME1'
1005 include 'COMMON.LAGRANGE'
1006 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1007 &,time01,zcopy(dimen3)
1010 if (nfgtasks.gt.1) then
1011 if (fg_rank.eq.0) then
1012 c The matching BROADCAST for fg processors is called in ERGASTULUM
1014 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
1015 time_Bcast=time_Bcast+MPI_Wtime()-time00
1016 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
1018 c write (2,*) "time00",time00
1019 c write (2,*) "Before Scatterv"
1021 c write (2,*) "Whole z (for FG master)"
1023 c write (2,*) i,z(i)
1025 c call MPI_Barrier(FG_COMM,IERROR)
1027 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1028 & MPI_DOUBLE_PRECISION,
1029 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1030 c write (2,*) "My chunk of z"
1031 do i=1,3*my_ng_count
1033 c write (2,*) i,z(i)
1035 c write (2,*) "After SCATTERV"
1037 c write (2,*) "MPI_Wtime",MPI_Wtime()
1038 time_scatter=time_scatter+MPI_Wtime()-time00
1040 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1042 c write (2,*) "time_scatter",time_scatter
1043 c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1052 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1053 c & Ginv(i,j),z((j-1)*3+k+1),
1054 c & Ginv(i,j)*z((j-1)*3+k+1)
1055 c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1056 temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
1060 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1061 c write (2,*) "Before REDUCE"
1063 c write (2,*) "z before reduce"
1065 c write (2,*) i,temp(i)
1068 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1069 & MPI_SUM,king,FG_COMM,IERR)
1070 time_reduce=time_reduce+MPI_Wtime()-time00
1071 c write (2,*) "After REDUCE"
1083 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1085 c & Ginv(i,j),z((j-1)*3+k+1),
1086 c & Ginv(i,j)*z((j-1)*3+k+1)
1087 d_a_tmp(ind)=d_a_tmp(ind)
1088 & +Ginv(j,i)*z((j-1)*3+k+1)
1089 c d_a_tmp(ind)=d_a_tmp(ind)
1090 c & +Ginv(i,j)*z((j-1)*3+k+1)
1095 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1102 c---------------------------------------------------------------------------
1104 SUBROUTINE ginv_mult_test(z,d_a_tmp)
1106 include 'DIMENSIONS'
1107 include 'COMMON.LAGRANGE'
1108 double precision z(dimen),d_a_tmp(dimen)
1109 double precision ztmp(dimen/3),dtmp(dimen/3)
1114 c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1123 ztmp(j)=z((j-1)*3+k+1)
1126 call alignx(16,ztmp(1))
1127 call alignx(16,dtmp(1))
1128 call alignx(16,Ginv(1,1))
1133 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1138 d_a_tmp(ind)=dtmp(i)
1144 c---------------------------------------------------------------------------
1145 SUBROUTINE fricmat_mult(z,d_a_tmp)
1146 include 'DIMENSIONS'
1152 include 'COMMON.LAGRANGE'
1153 include 'COMMON.IOUNITS'
1154 include 'COMMON.SETUP'
1155 include 'COMMON.TIME1'
1157 include 'COMMON.LANGEVIN'
1159 include 'COMMON.LANGEVIN.lang0'
1161 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1162 &,time01,zcopy(dimen3)
1164 if (nfgtasks.gt.1) then
1165 if (fg_rank.eq.0) then
1166 c The matching BROADCAST for fg processors is called in ERGASTULUM
1168 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1169 time_Bcast=time_Bcast+MPI_Wtime()-time00
1170 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1172 c call MPI_Barrier(FG_COMM,IERROR)
1174 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1175 & MPI_DOUBLE_PRECISION,
1176 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1177 c write (2,*) "My chunk of z"
1178 do i=1,3*my_ng_count
1180 c write (2,*) i,z(i)
1182 time_scatter=time_scatter+MPI_Wtime()-time00
1184 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1192 temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
1196 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1197 c write (2,*) "Before REDUCE"
1198 c write (2,*) "d_a_tmp before reduce"
1200 c write (2,*) i,temp(i)
1204 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1205 & MPI_SUM,king,FG_COMM,IERR)
1206 time_reduce=time_reduce+MPI_Wtime()-time00
1207 c write (2,*) "After REDUCE"
1219 d_a_tmp(ind)=d_a_tmp(ind)
1220 & -fricmat(j,i)*z((j-1)*3+k+1)
1225 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1230 c write (iout,*) "Vector d_a"
1232 c write (2,*) i,d_a_tmp(i)