2 !-----------------------------------------------------------------------------
9 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
14 !-----------------------------------------------------------------------------
16 !-----------------------------------------------------------------------------
18 !-------------------------------------------------------------------------
19 ! This subroutine contains the total lagrangain from which the accelerations
20 ! are obtained. For numerical gradient checking, the derivetive of the
21 ! lagrangian in the velocities and coordinates are calculated seperately
22 !-------------------------------------------------------------------------
23 ! implicit real*8 (a-h,o-z)
24 ! include 'DIMENSIONS'
27 use geometry_data, only: nres
28 use control_data !el, only: mucadyn,lmuca
31 real(kind=8) :: time00
33 ! include 'COMMON.VAR'
34 ! include 'COMMON.CHAIN'
35 ! include 'COMMON.DERIV'
36 ! include 'COMMON.GEO'
37 ! include 'COMMON.LOCAL'
38 ! include 'COMMON.INTERACT'
40 ! include 'COMMON.IOUNITS'
41 ! include 'COMMON.CONTROL'
42 ! include 'COMMON.MUCA'
43 ! include 'COMMON.TIME1'
44 integer ::i,j,ind,itime,k
45 real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres
46 real(kind=8) :: rs(dimen),xsolv(dimen)
48 real(kind=8) :: rscheck(dimen),rsold(dimen)
50 logical :: lprn = .false.
51 !el common /cipiszcze/ itime
59 write (iout,*) "Potential forces backbone"
61 write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
63 write (iout,*) "Potential forces sidechain"
65 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
66 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
72 if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
73 rs(ind)=-gcart(j,i)-gxcart(j,i)
77 rs(ind+1)=-gxcart(j,i)
85 write(iout,*) "RHS of the 5-diag equations system"
91 call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
93 write (iout,*) "Solution of the 5-diagonal equations system"
95 write (iout,'(i5,f10.5)') i,xsolv(i)
100 rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
102 rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
103 DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
106 rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
107 xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
108 xsolv(i+1)+DU2orig(i)*xsolv(i+2)
110 rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
111 DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
112 xsolv(dimen-1)+DU1orig(dimen-1)*&
114 rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
115 xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
118 write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
119 "ratio",rscheck(i)/rsold(i)
125 if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
130 d_a(j,i+nres)=xsolv(ind+1)
137 write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
140 write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
147 if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
149 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
153 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
154 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
162 write(iout,*) 'acceleration 3D FIVEDIAG'
163 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
165 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
168 write (iout,'(i3,3f10.5,3x,3f10.5)') &
169 i+nres,(d_a(j,i+nres),j=1,3)
179 write (iout,*) "Potential forces backbone"
182 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
183 i,(-gcart(j,i),j=1,3)
186 zapas(ind)=-gcart(j,i)
189 if (lprn) write (iout,*) "Potential forces sidechain"
191 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
192 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
193 i,(-gxcart(j,i),j=1,3)
196 zapas(ind)=-gxcart(j,i)
201 call ginv_mult(zapas,d_a_work)
210 d_a(j,i)=d_a_work(ind)
214 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
217 d_a(j,i+nres)=d_a_work(ind)
225 if(mucadyn.gt.0) call muca_update(potE)
226 factor=muca_factor(potE)*t_bath*Rb
228 !d print *,'lmuca ',factor,potE
230 d_a(j,0)=d_a(j,0)*factor
234 d_a(j,i)=d_a(j,i)*factor
239 d_a(j,i+nres)=d_a(j,i+nres)*factor
246 write(iout,*) 'acceleration 3D'
247 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
249 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
252 write (iout,'(i3,3f10.5,3x,3f10.5)') &
253 i+nres,(d_a(j,i+nres),j=1,3)
257 time_lagrangian=time_lagrangian+MPI_Wtime()-time00
260 end subroutine lagrangian
261 !-----------------------------------------------------------------------------
262 subroutine setup_MD_matrices
264 use geometry_data, only: nres,nside
268 use geometry, only:int_bounds
270 ! implicit real*8 (a-h,o-z)
271 ! include 'DIMENSIONS'
275 real(kind=8) :: time00
276 real(kind=8) ,allocatable, dimension(:) :: DDM,DDU1,DDU2
278 ! include 'COMMON.SETUP'
279 ! include 'COMMON.VAR'
280 ! include 'COMMON.CHAIN'
281 ! include 'COMMON.DERIV'
282 ! include 'COMMON.GEO'
283 ! include 'COMMON.LOCAL'
284 ! include 'COMMON.INTERACT'
285 ! include 'COMMON.MD'
287 ! include 'COMMON.LANGEVIN'
289 ! include 'COMMON.LANGEVIN.lang0'
291 ! include 'COMMON.IOUNITS'
292 ! include 'COMMON.TIME1'
293 logical :: lprn = .false.
296 real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
299 real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres
300 real(kind=8) :: relfeh,eps1,eps2
301 !el real(kind=8),dimension(:),allocatable :: Ghalf
302 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
303 !el real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
304 !el real(kind=8),dimension(:,:),allocatable :: Gcopy
305 real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
306 integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres
307 !el common /przechowalnia/ Gcopy,Ghalf
308 real(kind=8) :: coeff
309 integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
315 write (iout,*) "before FIVEDIAG"
317 write (iout,*) "ALLOCATE"
318 if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
319 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
322 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
323 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
325 ! Determine the number of degrees of freedom (dimen) and the number of
327 dimen=(nct-nnt+1)+nside
328 dimen1=(nct-nnt)+(nct-nnt+1)
330 write (iout,*) "nnt",nnt," nct",nct," nside",nside
334 if (iabs(itype(nnt)).eq.10) then
339 DM(1)=DM(1)+isc(iabs(itype(nnt)))
340 DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
345 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
347 if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
348 if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
351 DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
352 DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
358 if (iabs(itype(nct)).eq.10) then
359 DM(ind)=DM(ind)+msc(10)
362 DM(ind)=DM(ind)+isc(iabs(itype(nct)))
363 DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
371 if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
372 DU1(ind)=-isc(iabs(itype(i,1)))
383 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
384 write (iout,*) "i",i," itype",itype(i,1),ntyp1
385 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,1)).ne.ntyp1) then
398 write (iout,*) "nind",nind," dimen",dimen
399 if (nind.ne.dimen) then
400 write (iout,*) "ERROR, dimensions differ"
402 call MPI_Finalize(ierr)
406 write (iout,*) "The upper part of the five-diagonal inertia matrix"
408 if (i.lt.dimen-1) then
409 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
410 else if (i.eq.dimen-1) then
411 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
413 write (iout,'(i3,3f10.5)') i,DM(i)
417 call FDISYP (dimen, DM, DU1, DU2, MARK)
420 write (iout,*) "ERROR: the inertia matrix is not positive definite"
422 call MPI_Finalize(ierr)
425 else if (mark.eq.0) then
426 write (iout,*) "ERROR: the inertia matrix is singular"
428 call MPI_Finalize(ierr)
430 else if (mark.eq.1) then
431 write (iout,*) "The transformed inertia matrix"
433 if (i.lt.dimen-1) then
434 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
435 else if (i.eq.dimen-1) then
436 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
438 write (iout,'(i3,3f10.5)') i,DM(i)
442 ! Diagonalization of the pentadiagonal matrix
443 allocate(DDU1(2*nres))
444 allocate(DDU2(2*nres))
445 allocate(DDM(2*nres))
453 call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
456 "Eigenvalues of the five-diagonal inertia matrix"
458 write (iout,'(i5,f10.5)') i,geigen(i)
467 if (nfgtasks.gt.1) then
469 call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
470 time_Bcast=time_Bcast+MPI_Wtime()-time00
471 call int_bounds(dimen,igmult_start,igmult_end)
472 igmult_start=igmult_start-1
473 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
474 ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
475 my_ng_count=igmult_end-igmult_start
476 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
477 MPI_INTEGER,FG_COMM,IERROR)
478 write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
479 ' absolute rank',myrank,' igmult_start',igmult_start,&
480 ' igmult_end',igmult_end,' count',my_ng_count
481 write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
482 write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
492 ! write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
493 ! Zeroing out A and fricmat
499 ! Diagonal elements of the dC part of A and the respective friction coefficients
511 ! Off-diagonal elements of the dC part of A
518 ! Diagonal elements of the dX part of A and the respective friction coefficients
528 massvec(ii)=msc(iabs(iti),1)
529 if (iti.ne.10 .and. iti.ne.ntyp1) then
533 Gmat(ii1,ii1)=ISC(iabs(iti),1)
536 ! Off-diagonal elements of the dX part of A
550 write (iout,*) "Vector massvec"
552 write (iout,*) i,massvec(i)
554 write (iout,'(//a)') "A"
555 call matout(dimen,dimen1,nres2,nres2,A)
558 ! Calculate the G matrix (store in Gmat)
563 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
565 Gmat(k,i)=Gmat(k,i)+dtdi
570 write (iout,'(//a)') "Gmat"
571 call matout(dimen,dimen,nres2,nres2,Gmat)
580 ! Invert the G matrix
581 call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
583 write (iout,'(//a)') "Ginv"
584 call matout(dimen,dimen,nres2,nres2,Ginv)
591 if (itype(i,1).eq.10) then
599 Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
600 Bmat(i-nnt+1+(nct-nnt)+1,j)=1
603 Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
604 Bmat(i-nnt+1+(nct-nnt)+1,j)=1
608 write (iout,*) "j",j," dimen",dimen
610 write (iout,'(//a)') "Bmat"
611 call matout(dimen,dimen,nres2,nres2,Bmat)
618 Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
624 write (iout,'(//a)') "Gmat-transformed"
625 call matout(dimen,dimen,nres2,nres2,Gcopy)
628 if (nfgtasks.gt.1) then
629 myginv_ng_count=nres2*my_ng_count
630 call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
631 nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
632 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
633 nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
634 write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
635 write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
637 ! call MPI_Scatterv(ginv(1,1),nginv_counts(0),
638 ! & nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
639 ! & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
640 ! call MPI_Barrier(FG_COMM,IERR)
642 call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
643 nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
644 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
646 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
653 ! write (iout,*) "Master's chunk of ginv"
654 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
658 write (iout,*) "The G matrix is singular."
661 ! Compute G**(-1/2) and G**(1/2)
669 call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
672 write (iout,'(//a)') &
673 "Eigenvectors and eigenvalues of the G matrix"
674 call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
677 sqreig(i)=dsqrt(Geigen(i))
685 Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
686 Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
687 Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
692 write (iout,*) "Comparison of original and restored G"
695 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
696 Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
703 end subroutine setup_MD_matrices
704 !-----------------------------------------------------------------------------
705 subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
707 ! implicit real*8 (a-h,o-z)
708 ! include 'DIMENSIONS'
709 ! include 'COMMON.IOUNITS'
710 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
711 real(kind=8) :: A(LM2,LM3),B(LM2)
715 WRITE(IOUT,600) (I,I=KA,KB)
716 WRITE(IOUT,601) (B(I),I=KA,KB)
720 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
726 4 IF (KB.EQ.NC) RETURN
730 600 FORMAT (// 9H ROOT NO.,I4,9I11)
731 601 FORMAT (/5X,10(1PE11.4))
733 603 FORMAT (I5,10F11.5)
735 end subroutine EIGOUT
736 !-----------------------------------------------------------------------------
737 subroutine MATOUT(NC,NR,LM2,LM3,A)
739 ! implicit real*8 (a-h,o-z)
740 ! include 'DIMENSIONS'
741 ! include 'COMMON.IOUNITS'
742 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
743 real(kind=8) :: A(LM2,LM3)
747 WRITE(IOUT,600) (I,I=KA,KB)
751 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
757 4 IF (KB.EQ.NC) RETURN
761 600 FORMAT (//5x,9I11)
763 603 FORMAT (I5,10F11.3)
765 end subroutine MATOUT
766 !-----------------------------------------------------------------------------
767 subroutine MATOUT1(NC,NR,LM2,LM3,A)
769 ! implicit real*8 (a-h,o-z)
770 ! include 'DIMENSIONS'
771 ! include 'COMMON.IOUNITS'
772 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
773 real(kind=8) :: A(LM2,LM3)
777 WRITE(IOUT,600) (I,I=KA,KB)
781 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
787 4 IF (KB.EQ.NC) RETURN
791 600 FORMAT (//5x,7(3I5,2x))
793 603 FORMAT (I5,7(3F5.1,2x))
795 end subroutine MATOUT1
796 !-----------------------------------------------------------------------------
797 subroutine MATOUT2(NC,NR,LM2,LM3,A)
799 ! implicit real*8 (a-h,o-z)
800 ! include 'DIMENSIONS'
801 ! include 'COMMON.IOUNITS'
802 integer :: I,J,KA,KC,KB,N
803 integer :: LM2,LM3,NC,NR
804 real(kind=8) :: A(LM2,LM3)
808 WRITE(IOUT,600) (I,I=KA,KB)
812 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
818 4 IF (KB.EQ.NC) RETURN
822 600 FORMAT (//5x,4(3I9,2x))
824 603 FORMAT (I5,4(3F9.3,2x))
826 end subroutine MATOUT2
827 !-----------------------------------------------------------------------------
828 subroutine ginv_mult(z,d_a_tmp)
830 use geometry_data, only: nres
833 ! implicit real*8 (a-h,o-z)
834 ! include 'DIMENSIONS'
837 integer :: ierr,ierror
839 ! include 'COMMON.SETUP'
840 ! include 'COMMON.TIME1'
841 ! include 'COMMON.MD'
842 real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
843 real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres
844 real(kind=8) :: time00,time01
847 if (nfgtasks.gt.1) then
848 if (fg_rank.eq.0) then
849 ! The matching BROADCAST for fg processors is called in ERGASTULUM
851 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
852 time_Bcast=time_Bcast+MPI_Wtime()-time00
853 ! print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
855 ! write (2,*) "time00",time00
856 ! write (2,*) "Before Scatterv"
858 ! write (2,*) "Whole z (for FG master)"
862 ! call MPI_Barrier(FG_COMM,IERROR)
864 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
865 call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
866 MPI_DOUBLE_PRECISION,&
867 z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
868 ! write (2,*) "My chunk of z"
869 ! do i=1,3*my_ng_count
872 ! write (2,*) "After SCATTERV"
874 ! write (2,*) "MPI_Wtime",MPI_Wtime()
875 time_scatter=time_scatter+MPI_Wtime()-time00
877 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
879 ! write (2,*) "time_scatter",time_scatter
880 ! write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
889 ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
890 ! & Ginv(i,j),z((j-1)*3+k+1),
891 ! & Ginv(i,j)*z((j-1)*3+k+1)
892 ! temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
893 temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
897 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
898 ! write (2,*) "Before REDUCE"
900 ! write (2,*) "z before reduce"
902 ! write (2,*) i,temp(i)
905 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
906 MPI_SUM,king,FG_COMM,IERR)
907 time_reduce=time_reduce+MPI_Wtime()-time00
908 ! write (2,*) "After REDUCE"
920 ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
922 ! & Ginv(i,j),z((j-1)*3+k+1),
923 ! & Ginv(i,j)*z((j-1)*3+k+1)
924 d_a_tmp(ind)=d_a_tmp(ind) &
925 +Ginv(j,i)*z((j-1)*3+k+1)
926 ! d_a_tmp(ind)=d_a_tmp(ind)
927 ! & +Ginv(i,j)*z((j-1)*3+k+1)
932 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
938 end subroutine ginv_mult
939 !-----------------------------------------------------------------------------
941 subroutine ginv_mult_test(z,d_a_tmp)
943 ! include 'DIMENSIONS'
945 ! include 'COMMON.MD'
946 real(kind=8),dimension(dimen) :: z,d_a_tmp
947 real(kind=8),dimension(dimen/3) :: ztmp,dtmp
952 ! d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
961 ztmp(j)=z((j-1)*3+k+1)
964 call alignx(16,ztmp(1))
965 call alignx(16,dtmp(1))
966 call alignx(16,Ginv(1,1))
971 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
980 end subroutine ginv_mult_test
982 !-----------------------------------------------------------------------------
983 subroutine fricmat_mult(z,d_a_tmp)
985 use geometry_data, only: nres
988 ! include 'DIMENSIONS'
991 integer :: IERROR,ierr
993 ! include 'COMMON.MD'
994 ! include 'COMMON.IOUNITS'
995 ! include 'COMMON.SETUP'
996 ! include 'COMMON.TIME1'
998 ! include 'COMMON.LANGEVIN'
1000 ! include 'COMMON.LANGEVIN.lang0'
1002 real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1003 real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres
1004 real(kind=8) :: time00,time01
1005 integer :: i,j,k,ind,nres2
1007 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1010 if (nfgtasks.gt.1) then
1011 if (fg_rank.eq.0) then
1012 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1014 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1015 time_Bcast=time_Bcast+MPI_Wtime()-time00
1016 ! print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1018 ! call MPI_Barrier(FG_COMM,IERROR)
1020 call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1021 MPI_DOUBLE_PRECISION,&
1022 z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1023 ! write (2,*) "My chunk of z"
1024 ! do i=1,3*my_ng_count
1025 ! write (2,*) i,z(i)
1027 time_scatter=time_scatter+MPI_Wtime()-time00
1029 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1037 temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1041 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1042 ! write (2,*) "Before REDUCE"
1043 ! write (2,*) "d_a_tmp before reduce"
1045 ! write (2,*) i,temp(i)
1049 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1050 MPI_SUM,king,FG_COMM,IERR)
1051 time_reduce=time_reduce+MPI_Wtime()-time00
1052 ! write (2,*) "After REDUCE"
1064 d_a_tmp(ind)=d_a_tmp(ind) &
1065 -fricmat(j,i)*z((j-1)*3+k+1)
1070 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1075 ! write (iout,*) "Vector d_a"
1077 ! write (2,*) i,d_a_tmp(i)
1080 end subroutine fricmat_mult
1081 !-----------------------------------------------------------------------------
1082 !-----------------------------------------------------------------------------