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)
13 include 'COMMON.CHAIN'
14 include 'COMMON.DERIV'
16 include 'COMMON.LOCAL'
17 include 'COMMON.INTERACT'
19 include 'COMMON.IOUNITS'
20 include 'COMMON.CONTROL'
22 include 'COMMON.TIME1'
25 double precision zapas(MAXRES6),muca_factor
26 logical lprn /.false./
27 common /cipiszcze/ itime
37 write (iout,*) "Potential forces backbone"
40 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
41 & i,(-gcart(j,i),j=1,3)
44 zapas(ind)=-gcart(j,i)
47 if (lprn) write (iout,*) "Potential forces sidechain"
49 if (itype(i).ne.10) then
50 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
51 & i,(-gcart(j,i),j=1,3)
54 zapas(ind)=-gxcart(j,i)
59 call ginv_mult(zapas,d_a_work)
68 d_a(j,i)=d_a_work(ind)
72 if (itype(i).ne.10) then
75 d_a(j,i+nres)=d_a_work(ind)
82 if(mucadyn.gt.0) call muca_update(potE)
83 factor=muca_factor(potE)*t_bath*Rb
85 cd print *,'lmuca ',factor,potE
87 d_a(j,0)=d_a(j,0)*factor
91 d_a(j,i)=d_a(j,i)*factor
96 d_a(j,i+nres)=d_a(j,i+nres)*factor
103 write(iout,*) 'acceleration 3D'
104 write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
106 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
109 write (iout,'(i3,3f10.5,3x,3f10.5)')
110 & i+nres,(d_a(j,i+nres),j=1,3)
114 time_lagrangian=time_lagrangian+MPI_Wtime()-time00
118 c------------------------------------------------------------------
119 subroutine setup_MD_matrices
120 implicit real*8 (a-h,o-z)
126 include 'COMMON.SETUP'
127 include 'COMMON.CONTROL'
129 include 'COMMON.CHAIN'
130 include 'COMMON.DERIV'
132 include 'COMMON.LOCAL'
133 include 'COMMON.INTERACT'
136 include 'COMMON.LANGEVIN'
138 include 'COMMON.LANGEVIN.lang0'
140 include 'COMMON.IOUNITS'
141 include 'COMMON.TIME1'
143 logical lprn /.false./
145 double precision dtdi,massvec(maxres2),Gcopy(maxres2,maxres2),
146 & Ghalf(mmaxres2),sqreig(maxres2), invsqreig(maxres2), Gcopytmp,
147 & Gsqrptmp, Gsqrmtmp, Gvec2tmp,Gvectmp(maxres2,maxres2)
148 double precision work(8*maxres6)
149 integer iwork(maxres6)
150 common /przechowalnia/ Gcopy,Ghalf,invsqreig,Gvectmp
152 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
153 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
155 c Determine the number of degrees of freedom (dimen) and the number of
157 dimen=(nct-nnt+1)+nside
158 dimen1=(nct-nnt)+(nct-nnt+1)
161 if (nfgtasks.gt.1) then
163 call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
164 time_Bcast=time_Bcast+MPI_Wtime()-time00
165 call int_bounds(dimen,igmult_start,igmult_end)
166 igmult_start=igmult_start-1
167 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
168 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
169 my_ng_count=igmult_end-igmult_start
170 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
171 & MPI_INTEGER,FG_COMM,IERROR)
172 write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
173 & ' absolute rank',myrank,' igmult_start',igmult_start,
174 & ' igmult_end',igmult_end,' count',my_ng_count
175 write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
176 write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
186 c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
187 c Zeroing out A and fricmat
193 c Diagonal elements of the dC part of A and the respective friction coefficients
205 c Off-diagonal elements of the dC part of A
212 c Diagonal elements of the dX part of A and the respective friction coefficients
226 Gmat(ii1,ii1)=ISC(iti)
229 c Off-diagonal elements of the dX part of A
243 write (iout,*) "Vector massvec"
245 write (iout,*) i,massvec(i)
247 write (iout,'(//a)') "A"
248 call matout(dimen,dimen1,maxres2,maxres2,A)
251 c Calculate the G matrix (store in Gmat)
256 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
258 Gmat(k,i)=Gmat(k,i)+dtdi
263 write (iout,'(//a)') "Gmat"
264 call matout(dimen,dimen,maxres2,maxres2,Gmat)
273 c Invert the G matrix
274 call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
276 write (iout,'(//a)') "Ginv"
277 call matout(dimen,dimen,maxres2,maxres2,Ginv)
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 if (lprn .and. (me.eq.king .or. .not. out1file) ) then
287 write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
288 write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
291 c call MPI_Scatterv(ginv(1,1),nginv_counts(0),
292 c & nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
293 c & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
294 c call MPI_Barrier(FG_COMM,IERR)
296 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
297 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
298 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
300 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
307 c write (iout,*) "Master's chunk of ginv"
308 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
312 write (iout,*) "The G matrix is singular."
315 c Compute G**(-1/2) and G**(1/2)
323 call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
327 & "Eigenvectors and eigenvalues of the G matrix"
328 call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
332 sqreig(i)=dsqrt(Geigen(i))
333 invsqreig(i)=1.d0/sqreig(i)
337 Gvectmp(i,j)=Gvec(j,i)
347 c Gvec2tmp=Gvec(i,k)*Gvec(j,k)
348 Gvec2tmp=Gvec(k,i)*Gvec(k,j)
349 Gsqrptmp=Gsqrptmp+Gvec2tmp*sqreig(k)
350 Gsqrmtmp=Gsqrmtmp+Gvec2tmp*invsqreig(k)
351 Gcopytmp=Gcopytmp+Gvec2tmp*Geigen(k)
361 Gvec(i,j)=Gvectmp(j,i)
366 write (iout,*) "Comparison of original and restored G"
369 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
370 & Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
376 c-------------------------------------------------------------------------------
377 SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
378 implicit real*8 (a-h,o-z)
380 include 'COMMON.IOUNITS'
381 double precision A(LM2,LM3),B(LM2)
385 WRITE(IOUT,600) (I,I=KA,KB)
386 WRITE(IOUT,601) (B(I),I=KA,KB)
390 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
396 4 IF (KB.EQ.NC) RETURN
400 600 FORMAT (// 9H ROOT NO.,I4,9I11)
401 601 FORMAT (/5X,10(1PE11.4))
403 603 FORMAT (I5,10F11.5)
406 c-------------------------------------------------------------------------------
407 SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
408 implicit real*8 (a-h,o-z)
410 include 'COMMON.IOUNITS'
411 double precision A(LM2,LM3)
415 WRITE(IOUT,600) (I,I=KA,KB)
419 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
425 4 IF (KB.EQ.NC) RETURN
429 600 FORMAT (//5x,9I11)
431 603 FORMAT (I5,10F11.3)
434 c-------------------------------------------------------------------------------
435 SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
436 implicit real*8 (a-h,o-z)
438 include 'COMMON.IOUNITS'
439 double precision A(LM2,LM3)
443 WRITE(IOUT,600) (I,I=KA,KB)
447 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
453 4 IF (KB.EQ.NC) RETURN
457 600 FORMAT (//5x,7(3I5,2x))
459 603 FORMAT (I5,7(3F5.1,2x))
462 c-------------------------------------------------------------------------------
463 SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
464 implicit real*8 (a-h,o-z)
466 include 'COMMON.IOUNITS'
467 double precision A(LM2,LM3)
471 WRITE(IOUT,600) (I,I=KA,KB)
475 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
481 4 IF (KB.EQ.NC) RETURN
485 600 FORMAT (//5x,4(3I9,2x))
487 603 FORMAT (I5,4(3F9.3,2x))
490 c---------------------------------------------------------------------------
491 SUBROUTINE ginv_mult(z,d_a_tmp)
492 implicit real*8 (a-h,o-z)
498 include 'COMMON.SETUP'
499 include 'COMMON.TIME1'
501 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
504 if (nfgtasks.gt.1) then
505 if (fg_rank.eq.0) then
506 c The matching BROADCAST for fg processors is called in ERGASTULUM
508 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
509 time_Bcast=time_Bcast+MPI_Wtime()-time00
510 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
512 c write (2,*) "time00",time00
513 c write (2,*) "Before Scatterv"
515 c write (2,*) "Whole z (for FG master)"
519 c call MPI_Barrier(FG_COMM,IERROR)
521 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
522 & MPI_DOUBLE_PRECISION,
523 & z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
524 c write (2,*) "My chunk of z"
525 c do i=1,3*my_ng_count
528 c write (2,*) "After SCATTERV"
530 c write (2,*) "MPI_Wtime",MPI_Wtime()
531 time_scatter=time_scatter+MPI_Wtime()-time00
533 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
535 c write (2,*) "time_scatter",time_scatter
536 c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
545 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
546 c & Ginv(i,j),z((j-1)*3+k+1),
547 c & Ginv(i,j)*z((j-1)*3+k+1)
548 c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
549 temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
553 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
554 c write (2,*) "Before REDUCE"
556 c write (2,*) "z before reduce"
558 c write (2,*) i,temp(i)
561 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
562 & MPI_SUM,king,FG_COMM,IERR)
563 time_reduce=time_reduce+MPI_Wtime()-time00
564 c write (2,*) "After REDUCE"
576 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
578 c & Ginv(i,j),z((j-1)*3+k+1),
579 c & Ginv(i,j)*z((j-1)*3+k+1)
580 d_a_tmp(ind)=d_a_tmp(ind)
581 & +Ginv(j,i)*z((j-1)*3+k+1)
582 c d_a_tmp(ind)=d_a_tmp(ind)
583 c & +Ginv(i,j)*z((j-1)*3+k+1)
588 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
595 c---------------------------------------------------------------------------
597 SUBROUTINE ginv_mult_test(z,d_a_tmp)
600 c include 'COMMON.MD'
601 double precision z(dimen),d_a_tmp(dimen)
602 double precision ztmp(dimen/3),dtmp(dimen/3)
607 c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
616 ztmp(j)=z((j-1)*3+k+1)
619 call alignx(16,ztmp(1))
620 call alignx(16,dtmp(1))
621 call alignx(16,Ginv(1,1))
626 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
637 c---------------------------------------------------------------------------
638 SUBROUTINE fricmat_mult(z,d_a_tmp)
645 include 'COMMON.IOUNITS'
646 include 'COMMON.SETUP'
647 include 'COMMON.TIME1'
649 include 'COMMON.LANGEVIN'
651 include 'COMMON.LANGEVIN.lang0'
653 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
656 if (nfgtasks.gt.1) then
657 if (fg_rank.eq.0) then
658 c The matching BROADCAST for fg processors is called in ERGASTULUM
660 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
661 time_Bcast=time_Bcast+MPI_Wtime()-time00
662 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
664 c call MPI_Barrier(FG_COMM,IERROR)
666 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
667 & MPI_DOUBLE_PRECISION,
668 & z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
669 c write (2,*) "My chunk of z"
670 c do i=1,3*my_ng_count
673 time_scatter=time_scatter+MPI_Wtime()-time00
675 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
683 temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
687 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
688 c write (2,*) "Before REDUCE"
689 c write (2,*) "d_a_tmp before reduce"
691 c write (2,*) i,temp(i)
695 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
696 & MPI_SUM,king,FG_COMM,IERR)
697 time_reduce=time_reduce+MPI_Wtime()-time00
698 c write (2,*) "After REDUCE"
710 d_a_tmp(ind)=d_a_tmp(ind)
711 & -fricmat(j,i)*z((j-1)*3+k+1)
716 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
721 c write (iout,*) "Vector d_a"
723 c write (2,*) i,d_a_tmp(i)