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,mnum
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"
66 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.lt.4)&
67 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
74 if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.ge.4)then
75 rs(ind)=-gcart(j,i)-gxcart(j,i)
79 rs(ind+1)=-gxcart(j,i)
87 write(iout,*) "RHS of the 5-diag equations system"
93 call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
95 write (iout,*) "Solution of the 5-diagonal equations system"
97 write (iout,'(i5,f10.5)') i,xsolv(i)
102 rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
104 rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
105 DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
108 rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
109 xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
110 xsolv(i+1)+DU2orig(i)*xsolv(i+2)
112 rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
113 DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
114 xsolv(dimen-1)+DU1orig(dimen-1)*&
116 rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
117 xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
120 write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
121 "ratio",rscheck(i)/rsold(i)
128 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or. mnum.ge.4)then
133 d_a(j,i+nres)=xsolv(ind+1)
140 write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
143 write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
151 ! if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
152 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or. mnum.ge.4)then
154 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
158 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
159 d_a(j,i)=d_a(j,i+1)-d_a(j,i)
167 write(iout,*) 'acceleration 3D FIVEDIAG'
168 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
170 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
173 write (iout,'(i3,3f10.5,3x,3f10.5)') &
174 i+nres,(d_a(j,i+nres),j=1,3)
184 write (iout,*) "Potential forces backbone"
187 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
188 i,(-gcart(j,i),j=1,3)
191 zapas(ind)=-gcart(j,i)
194 if (lprn) write (iout,*) "Potential forces sidechain"
197 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
199 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
200 i,(-gxcart(j,i),j=1,3)
203 zapas(ind)=-gxcart(j,i)
208 call ginv_mult(zapas,d_a_work)
217 d_a(j,i)=d_a_work(ind)
222 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
223 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
227 d_a(j,i+nres)=d_a_work(ind)
235 if(mucadyn.gt.0) call muca_update(potE)
236 factor=muca_factor(potE)*t_bath*Rb
238 !d print *,'lmuca ',factor,potE
240 d_a(j,0)=d_a(j,0)*factor
244 d_a(j,i)=d_a(j,i)*factor
249 d_a(j,i+nres)=d_a(j,i+nres)*factor
256 write(iout,*) 'acceleration 3D'
257 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
259 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
262 write (iout,'(i3,3f10.5,3x,3f10.5)') &
263 i+nres,(d_a(j,i+nres),j=1,3)
267 time_lagrangian=time_lagrangian+MPI_Wtime()-time00
270 end subroutine lagrangian
271 !-----------------------------------------------------------------------------
272 subroutine setup_MD_matrices
274 use geometry_data, only: nres,nside
278 use geometry, only:int_bounds
280 ! implicit real*8 (a-h,o-z)
281 ! include 'DIMENSIONS'
285 real(kind=8) :: time00
286 real(kind=8) ,allocatable, dimension(:) :: DDM,DDU1,DDU2
288 ! include 'COMMON.SETUP'
289 ! include 'COMMON.VAR'
290 ! include 'COMMON.CHAIN'
291 ! include 'COMMON.DERIV'
292 ! include 'COMMON.GEO'
293 ! include 'COMMON.LOCAL'
294 ! include 'COMMON.INTERACT'
295 ! include 'COMMON.MD'
297 ! include 'COMMON.LANGEVIN'
299 ! include 'COMMON.LANGEVIN.lang0'
301 ! include 'COMMON.IOUNITS'
302 ! include 'COMMON.TIME1'
303 logical :: lprn = .false.
306 real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
309 real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres
310 real(kind=8) :: relfeh,eps1,eps2
311 !el real(kind=8),dimension(:),allocatable :: Ghalf
312 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
313 !el real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
314 !el real(kind=8),dimension(:,:),allocatable :: Gcopy
315 real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
316 integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres
317 !el common /przechowalnia/ Gcopy,Ghalf
318 real(kind=8) :: coeff,mscab
319 integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
322 print *,"just entered"
326 write (iout,*) "before FIVEDIAG"
328 write (iout,*) "ALLOCATE"
330 if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
331 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
332 if(.not.allocated(Bmat)) allocate(Bmat(nres2,nres2))
333 if(.not.allocated(matmult)) allocate(matmult(nres2,nres2))
336 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
337 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
339 ! Determine the number of degrees of freedom (dimen) and the number of
341 dimen=(nct-nnt+1)+nside
342 dimen1=(nct-nnt)+(nct-nnt+1)
344 write (iout,*) "nnt",nnt," nct",nct," nside",nside
349 if (iabs(itype(nnt)).eq.10) then
354 DM(1)=DM(1)+isc(iabs(itype(nnt)))
355 DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
361 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
363 if (iabs(itype(i,1)).eq.10 .or. &
364 iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then
365 if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
368 DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
369 DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
375 if (iabs(itype(nct)).eq.10) then
376 DM(ind)=DM(ind)+msc(10)
379 DM(ind)=DM(ind)+isc(iabs(itype(nct)))
380 DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
389 if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
391 DU1(ind)=-isc(iabs(itype(i,1)))
403 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
404 write (iout,*) "i",i," itype",itype(i,1),ntyp1
405 if (iabs(itype(i,1)).ne.10 .and. &
406 iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
420 if (iabs(itype(nnt,1)).eq.10) then
421 DM(1)=DM(1)+msc(10,1)
425 DM(1)=DM(1)+isc(iabs(itype(nnt,mnum)),mnum)
426 DM(2)=msc(iabs(itype(nnt,mnum)),mnum)+isc(iabs(itype(nnt,mnum)),mnum)
431 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
433 DM(ind)=2*ip4+mp(1)/2
434 if (iabs(itype(i,1)).eq.10 .or. &
435 iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then
436 if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1)
437 if (mnum.eq.5) DM(ind)=DM(ind)+msc(itype(i,mnum),mnum)
440 DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum)
441 DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum)
447 if (iabs(itype(nct,1)).eq.10) then
448 DM(ind)=DM(ind)+msc(10,1)
452 DM(ind)=DM(ind)+isc(iabs(itype(nct,mnum)),mnum)
453 DM(ind+1)=msc(iabs(itype(nct,mnum)),mnum)+isc(iabs(itype(nct,mnum)),mnum)
462 ! if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1 &
463 ! .and.mnum.eq.5) then
464 if (iabs(itype(i,1)).ne.10 .and. &
465 iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
466 DU1(ind)=-isc(iabs(itype(i,1)),1)
471 DU1(ind)=msc(itype(i,mnum),5))
481 ! if (iabs(itype(i,1)).eq.ntyp1) cycle
482 write (iout,*) "i",i," itype",itype(i,1),ntyp1
483 if (iabs(itype(i,1)).ne.10 .and. &
484 iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
498 write (iout,*) "nind",nind," dimen",dimen
499 if (nind.ne.dimen) then
500 write (iout,*) "ERROR, dimensions differ"
502 call MPI_Finalize(ierr)
506 write (iout,*) "The upper part of the five-diagonal inertia matrix"
508 if (i.lt.dimen-1) then
509 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
510 else if (i.eq.dimen-1) then
511 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
513 write (iout,'(i3,3f10.5)') i,DM(i)
517 call FDISYP (dimen, DM, DU1, DU2, MARK)
520 write (iout,*) "ERROR: the inertia matrix is not positive definite"
522 call MPI_Finalize(ierr)
525 else if (mark.eq.0) then
526 write (iout,*) "ERROR: the inertia matrix is singular"
528 call MPI_Finalize(ierr)
530 else if (mark.eq.1) then
531 write (iout,*) "The transformed inertia matrix"
533 if (i.lt.dimen-1) then
534 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
535 else if (i.eq.dimen-1) then
536 write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
538 write (iout,'(i3,3f10.5)') i,DM(i)
542 ! Diagonalization of the pentadiagonal matrix
543 allocate(DDU1(2*nres))
544 allocate(DDU2(2*nres))
545 allocate(DDM(2*nres))
553 call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
556 "Eigenvalues of the five-diagonal inertia matrix"
558 write (iout,'(i5,f10.5)') i,geigen(i)
561 if (allocated(DDU1)) deallocate(DDU1)
562 if (allocated(DDU2)) deallocate(DDU2)
563 if (allocated(DDM)) deallocate(DDM)
567 if (nfgtasks.gt.1) then
569 call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
570 time_Bcast=time_Bcast+MPI_Wtime()-time00
571 call int_bounds(dimen,igmult_start,igmult_end)
572 igmult_start=igmult_start-1
573 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
574 ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
575 my_ng_count=igmult_end-igmult_start
576 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
577 MPI_INTEGER,FG_COMM,IERROR)
578 write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
579 ' absolute rank',myrank,' igmult_start',igmult_start,&
580 ' igmult_end',igmult_end,' count',my_ng_count
581 write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
582 write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
592 ! write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
593 ! Zeroing out A and fricmat
599 ! Diagonal elements of the dC part of A and the respective friction coefficients
602 ! print *,"TUTUTUT?!",nnt,nct-1
607 coeff=0.25d0*IP(mnum)
608 print *,"TU",i,itype(i,mnum),mnum
609 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
610 print *,"TU2",i,coeff,mp(mnum)
611 massvec(ind1)=mp(mnum)
613 ! print *,"i",mp(mnum)
617 ! Off-diagonal elements of the dC part of A
624 ! Diagonal elements of the dX part of A and the respective friction coefficients
630 msc(ntyp1_molec(i),i)=1.0d0
632 msc(ntyp1_molec(4),4)=1.0d0
633 ! print *,"TU3",ntyp1_molec(4)-1
634 msc(ntyp1_molec(4)-1,4)=1.0d0
640 ! if (mnum.eq.5) then
643 mscab=msc(iabs(iti),mnum)
647 if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
651 Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
654 ! Off-diagonal elements of the dX part of A
669 write (iout,*) "Vector massvec"
671 write (iout,*) i,massvec(i)
673 write (iout,'(//a)') "A"
674 call matout(dimen,dimen1,nres2,nres2,A)
677 ! Calculate the G matrix (store in Gmat)
682 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
684 Gmat(k,i)=Gmat(k,i)+dtdi
689 write (iout,'(//a)') "Gmat"
690 call matout(dimen,dimen,nres2,nres2,Gmat)
699 ! Invert the G matrix
700 call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
702 write (iout,'(//a)') "Ginv"
703 call matout(dimen,dimen,nres2,nres2,Ginv)
711 if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
712 .or.(mnum.ge.4)) then
720 Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
721 Bmat(i-nnt+1+(nct-nnt)+1,j)=1
724 Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
725 Bmat(i-nnt+1+(nct-nnt)+1,j)=1
729 write (iout,*) "j",j," dimen",dimen
731 write (iout,'(//a)') "Bmat"
732 call matout(dimen,dimen,nres2,nres2,Bmat)
735 write(iout,*) "before Gcopy",dimen,nres*2
739 ! write (iout,*) "myij",i,j
742 ! write(iout,*) "i,j,k,l",i,j,k,l
743 Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
750 write (iout,'(//a)') "Gmat-transformed"
751 call matout(dimen,dimen,nres2,nres2,Gcopy)
754 if (nfgtasks.gt.1) then
755 myginv_ng_count=nres2*my_ng_count
756 call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
757 nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
758 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
759 nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
760 write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
761 write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
763 ! call MPI_Scatterv(ginv(1,1),nginv_counts(0),
764 ! & nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
765 ! & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
766 ! call MPI_Barrier(FG_COMM,IERR)
768 call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
769 nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
770 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
772 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
779 ! write (iout,*) "Master's chunk of ginv"
780 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
784 write (iout,*) "The G matrix is singular."
785 write (iout,'(//a)') "Gmat-transformed"
786 call matout(dimen,dimen,nres2,nres2,Gcopy)
789 ! Compute G**(-1/2) and G**(1/2)
797 call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
800 write (iout,'(//a)') &
801 "Eigenvectors and eigenvalues of the G matrix"
802 call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
805 sqreig(i)=dsqrt(Geigen(i))
813 Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
814 Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
815 Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
820 write (iout,*) "Comparison of original and restored G"
823 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
824 Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
828 if (allocated(Gcopy)) deallocate(Gcopy)
830 !write(iout,*) "end setup_MD_matr"
832 end subroutine setup_MD_matrices
833 !-----------------------------------------------------------------------------
834 subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
836 ! implicit real*8 (a-h,o-z)
837 ! include 'DIMENSIONS'
838 ! include 'COMMON.IOUNITS'
839 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
840 real(kind=8) :: A(LM2,LM3),B(LM2)
844 WRITE(IOUT,600) (I,I=KA,KB)
845 WRITE(IOUT,601) (B(I),I=KA,KB)
849 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
855 4 IF (KB.EQ.NC) RETURN
859 600 FORMAT (// 9H ROOT NO.,I4,9I11)
860 601 FORMAT (/5X,10(1PE11.4))
862 603 FORMAT (I5,10F11.5)
864 end subroutine EIGOUT
865 !-----------------------------------------------------------------------------
866 subroutine MATOUT(NC,NR,LM2,LM3,A)
868 ! implicit real*8 (a-h,o-z)
869 ! include 'DIMENSIONS'
870 ! include 'COMMON.IOUNITS'
871 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
872 real(kind=8) :: A(LM2,LM3)
876 WRITE(IOUT,600) (I,I=KA,KB)
880 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
886 4 IF (KB.EQ.NC) RETURN
890 600 FORMAT (//5x,9I11)
892 603 FORMAT (I5,10F11.3)
894 end subroutine MATOUT
895 !-----------------------------------------------------------------------------
896 subroutine MATOUT1(NC,NR,LM2,LM3,A)
898 ! implicit real*8 (a-h,o-z)
899 ! include 'DIMENSIONS'
900 ! include 'COMMON.IOUNITS'
901 integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
902 real(kind=8) :: A(LM2,LM3)
906 WRITE(IOUT,600) (I,I=KA,KB)
910 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
916 4 IF (KB.EQ.NC) RETURN
920 600 FORMAT (//5x,7(3I5,2x))
922 603 FORMAT (I5,7(3F5.1,2x))
924 end subroutine MATOUT1
925 !-----------------------------------------------------------------------------
926 subroutine MATOUT2(NC,NR,LM2,LM3,A)
928 ! implicit real*8 (a-h,o-z)
929 ! include 'DIMENSIONS'
930 ! include 'COMMON.IOUNITS'
931 integer :: I,J,KA,KC,KB,N
932 integer :: LM2,LM3,NC,NR
933 real(kind=8) :: A(LM2,LM3)
937 WRITE(IOUT,600) (I,I=KA,KB)
941 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
947 4 IF (KB.EQ.NC) RETURN
951 600 FORMAT (//5x,4(3I9,2x))
953 603 FORMAT (I5,4(3F9.3,2x))
955 end subroutine MATOUT2
956 !-----------------------------------------------------------------------------
957 subroutine ginv_mult(z,d_a_tmp)
959 use geometry_data, only: nres
962 ! implicit real*8 (a-h,o-z)
963 ! include 'DIMENSIONS'
966 integer :: ierr,ierror
968 ! include 'COMMON.SETUP'
969 ! include 'COMMON.TIME1'
970 ! include 'COMMON.MD'
971 real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
972 real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres
973 real(kind=8) :: time00,time01
976 if (nfgtasks.gt.1) then
977 if (fg_rank.eq.0) then
978 ! The matching BROADCAST for fg processors is called in ERGASTULUM
980 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
981 time_Bcast=time_Bcast+MPI_Wtime()-time00
982 ! print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
984 ! write (2,*) "time00",time00
985 ! write (2,*) "Before Scatterv"
987 ! write (2,*) "Whole z (for FG master)"
991 ! call MPI_Barrier(FG_COMM,IERROR)
993 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
994 call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
995 MPI_DOUBLE_PRECISION,&
996 z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
997 ! write (2,*) "My chunk of z"
998 ! do i=1,3*my_ng_count
1001 ! write (2,*) "After SCATTERV"
1003 ! write (2,*) "MPI_Wtime",MPI_Wtime()
1004 time_scatter=time_scatter+MPI_Wtime()-time00
1006 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1008 ! write (2,*) "time_scatter",time_scatter
1009 ! write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1018 ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1019 ! & Ginv(i,j),z((j-1)*3+k+1),
1020 ! & Ginv(i,j)*z((j-1)*3+k+1)
1021 ! temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1022 temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
1026 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1027 ! write (2,*) "Before REDUCE"
1029 ! write (2,*) "z before reduce"
1031 ! write (2,*) i,temp(i)
1034 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1035 MPI_SUM,king,FG_COMM,IERR)
1036 time_reduce=time_reduce+MPI_Wtime()-time00
1037 ! write (2,*) "After REDUCE"
1049 ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1051 ! & Ginv(i,j),z((j-1)*3+k+1),
1052 ! & Ginv(i,j)*z((j-1)*3+k+1)
1053 d_a_tmp(ind)=d_a_tmp(ind) &
1054 +Ginv(j,i)*z((j-1)*3+k+1)
1055 ! d_a_tmp(ind)=d_a_tmp(ind)
1056 ! & +Ginv(i,j)*z((j-1)*3+k+1)
1061 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1067 end subroutine ginv_mult
1068 !-----------------------------------------------------------------------------
1070 subroutine ginv_mult_test(z,d_a_tmp)
1072 ! include 'DIMENSIONS'
1073 !el integer :: dimen
1074 ! include 'COMMON.MD'
1075 real(kind=8),dimension(dimen) :: z,d_a_tmp
1076 real(kind=8),dimension(dimen/3) :: ztmp,dtmp
1077 integer :: i,j,k,ind
1081 ! d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1090 ztmp(j)=z((j-1)*3+k+1)
1093 call alignx(16,ztmp(1))
1094 call alignx(16,dtmp(1))
1095 call alignx(16,Ginv(1,1))
1100 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1105 d_a_tmp(ind)=dtmp(i)
1109 end subroutine ginv_mult_test
1111 !-----------------------------------------------------------------------------
1112 subroutine fricmat_mult(z,d_a_tmp)
1114 use geometry_data, only: nres
1117 ! include 'DIMENSIONS'
1120 integer :: IERROR,ierr
1122 ! include 'COMMON.MD'
1123 ! include 'COMMON.IOUNITS'
1124 ! include 'COMMON.SETUP'
1125 ! include 'COMMON.TIME1'
1127 ! include 'COMMON.LANGEVIN'
1129 ! include 'COMMON.LANGEVIN.lang0'
1131 real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
1132 real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres
1133 real(kind=8) :: time00,time01
1134 integer :: i,j,k,ind,nres2
1136 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
1139 if (nfgtasks.gt.1) then
1140 if (fg_rank.eq.0) then
1141 ! The matching BROADCAST for fg processors is called in ERGASTULUM
1143 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1144 time_Bcast=time_Bcast+MPI_Wtime()-time00
1145 ! print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1147 ! call MPI_Barrier(FG_COMM,IERROR)
1149 call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
1150 MPI_DOUBLE_PRECISION,&
1151 z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1152 ! write (2,*) "My chunk of z"
1153 ! do i=1,3*my_ng_count
1154 ! write (2,*) i,z(i)
1156 time_scatter=time_scatter+MPI_Wtime()-time00
1158 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1166 temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
1170 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1171 ! write (2,*) "Before REDUCE"
1172 ! write (2,*) "d_a_tmp before reduce"
1174 ! write (2,*) i,temp(i)
1178 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
1179 MPI_SUM,king,FG_COMM,IERR)
1180 time_reduce=time_reduce+MPI_Wtime()-time00
1181 ! write (2,*) "After REDUCE"
1193 d_a_tmp(ind)=d_a_tmp(ind) &
1194 -fricmat(j,i)*z((j-1)*3+k+1)
1199 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1204 ! write (iout,*) "Vector d_a"
1206 ! write (2,*) i,d_a_tmp(i)
1209 end subroutine fricmat_mult
1210 !-----------------------------------------------------------------------------
1211 !-----------------------------------------------------------------------------