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 write (iout,*) "The number of degrees of freedom ",dimen3
188 c Zeroing out A and fricmat
194 c Diagonal elements of the dC part of A and the respective friction coefficients
206 c Off-diagonal elements of the dC part of A
213 c Diagonal elements of the dX part of A and the respective friction coefficients
227 Gmat(ii1,ii1)=ISC(iti)
230 c Off-diagonal elements of the dX part of A
244 write (iout,*) "Vector massvec"
246 write (iout,*) i,massvec(i)
248 write (iout,'(//a)') "A"
249 call matout(dimen,dimen1,maxres2,maxres2,A)
252 c Calculate the G matrix (store in Gmat)
257 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
259 Gmat(k,i)=Gmat(k,i)+dtdi
264 write (iout,'(//a)') "Gmat"
265 call matout(dimen,dimen,maxres2,maxres2,Gmat)
274 c Invert the G matrix
275 call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
277 write (iout,'(//a)') "Ginv"
278 call matout(dimen,dimen,maxres2,maxres2,Ginv)
281 if (nfgtasks.gt.1) then
282 myginv_ng_count=maxres2*my_ng_count
283 call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
284 & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
285 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
286 & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
287 if (lprn .and. (me.eq.king .or. .not. out1file) ) then
288 write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
289 write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
292 c call MPI_Scatterv(ginv(1,1),nginv_counts(0),
293 c & nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
294 c & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
295 c call MPI_Barrier(FG_COMM,IERR)
297 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
298 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
299 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
301 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
308 c write (iout,*) "Master's chunk of ginv"
309 c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
313 write (iout,*) "The G matrix is singular."
316 c Compute G**(-1/2) and G**(1/2)
324 call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
328 & "Eigenvectors and eigenvalues of the G matrix"
329 call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
333 sqreig(i)=dsqrt(Geigen(i))
334 invsqreig(i)=1.d0/sqreig(i)
338 Gvectmp(i,j)=Gvec(j,i)
348 c Gvec2tmp=Gvec(i,k)*Gvec(j,k)
349 Gvec2tmp=Gvec(k,i)*Gvec(k,j)
350 Gsqrptmp=Gsqrptmp+Gvec2tmp*sqreig(k)
351 Gsqrmtmp=Gsqrmtmp+Gvec2tmp*invsqreig(k)
352 Gcopytmp=Gcopytmp+Gvec2tmp*Geigen(k)
362 Gvec(i,j)=Gvectmp(j,i)
367 write (iout,*) "Comparison of original and restored G"
370 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
371 & Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
377 c-------------------------------------------------------------------------------
378 SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
379 implicit real*8 (a-h,o-z)
381 include 'COMMON.IOUNITS'
382 double precision A(LM2,LM3),B(LM2)
386 WRITE(IOUT,600) (I,I=KA,KB)
387 WRITE(IOUT,601) (B(I),I=KA,KB)
391 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
397 4 IF (KB.EQ.NC) RETURN
401 600 FORMAT (// 9H ROOT NO.,I4,9I11)
402 601 FORMAT (/5X,10(1PE11.4))
404 603 FORMAT (I5,10F11.5)
407 c-------------------------------------------------------------------------------
408 SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
409 implicit real*8 (a-h,o-z)
411 include 'COMMON.IOUNITS'
412 double precision A(LM2,LM3)
416 WRITE(IOUT,600) (I,I=KA,KB)
420 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
426 4 IF (KB.EQ.NC) RETURN
430 600 FORMAT (//5x,9I11)
432 603 FORMAT (I5,10F11.3)
435 c-------------------------------------------------------------------------------
436 SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
437 implicit real*8 (a-h,o-z)
439 include 'COMMON.IOUNITS'
440 double precision A(LM2,LM3)
444 WRITE(IOUT,600) (I,I=KA,KB)
448 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
454 4 IF (KB.EQ.NC) RETURN
458 600 FORMAT (//5x,7(3I5,2x))
460 603 FORMAT (I5,7(3F5.1,2x))
463 c-------------------------------------------------------------------------------
464 SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
465 implicit real*8 (a-h,o-z)
467 include 'COMMON.IOUNITS'
468 double precision A(LM2,LM3)
472 WRITE(IOUT,600) (I,I=KA,KB)
476 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
482 4 IF (KB.EQ.NC) RETURN
486 600 FORMAT (//5x,4(3I9,2x))
488 603 FORMAT (I5,4(3F9.3,2x))
491 c---------------------------------------------------------------------------
492 SUBROUTINE ginv_mult(z,d_a_tmp)
493 implicit real*8 (a-h,o-z)
499 include 'COMMON.SETUP'
500 include 'COMMON.TIME1'
502 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
503 &,time01,zcopy(dimen3)
505 if (nfgtasks.gt.1) then
506 if (fg_rank.eq.0) then
507 c The matching BROADCAST for fg processors is called in ERGASTULUM
509 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
510 time_Bcast=time_Bcast+MPI_Wtime()-time00
511 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
513 c write (2,*) "time00",time00
514 c write (2,*) "Before Scatterv"
516 c write (2,*) "Whole z (for FG master)"
520 c call MPI_Barrier(FG_COMM,IERROR)
522 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
523 & MPI_DOUBLE_PRECISION,
524 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
528 c write (2,*) "My chunk of z"
529 c do i=1,3*my_ng_count
532 c write (2,*) "After SCATTERV"
534 c write (2,*) "MPI_Wtime",MPI_Wtime()
535 time_scatter=time_scatter+MPI_Wtime()-time00
537 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
539 c write (2,*) "time_scatter",time_scatter
540 c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
549 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
550 c & Ginv(i,j),z((j-1)*3+k+1),
551 c & Ginv(i,j)*z((j-1)*3+k+1)
552 c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
553 temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
557 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
558 c write (2,*) "Before REDUCE"
560 c write (2,*) "z before reduce"
562 c write (2,*) i,temp(i)
565 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
566 & MPI_SUM,king,FG_COMM,IERR)
567 time_reduce=time_reduce+MPI_Wtime()-time00
568 c write (2,*) "After REDUCE"
580 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
582 c & Ginv(i,j),z((j-1)*3+k+1),
583 c & Ginv(i,j)*z((j-1)*3+k+1)
584 d_a_tmp(ind)=d_a_tmp(ind)
585 & +Ginv(j,i)*z((j-1)*3+k+1)
586 c d_a_tmp(ind)=d_a_tmp(ind)
587 c & +Ginv(i,j)*z((j-1)*3+k+1)
592 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
599 c---------------------------------------------------------------------------
601 SUBROUTINE ginv_mult_test(z,d_a_tmp)
604 c include 'COMMON.MD'
605 double precision z(dimen),d_a_tmp(dimen)
606 double precision ztmp(dimen/3),dtmp(dimen/3)
611 c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
620 ztmp(j)=z((j-1)*3+k+1)
623 call alignx(16,ztmp(1))
624 call alignx(16,dtmp(1))
625 call alignx(16,Ginv(1,1))
630 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
641 c---------------------------------------------------------------------------
642 SUBROUTINE fricmat_mult(z,d_a_tmp)
649 include 'COMMON.IOUNITS'
650 include 'COMMON.SETUP'
651 include 'COMMON.TIME1'
653 include 'COMMON.LANGEVIN'
655 include 'COMMON.LANGEVIN.lang0'
657 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
658 &,time01,zcopy(dimen3)
660 if (nfgtasks.gt.1) then
661 if (fg_rank.eq.0) then
662 c The matching BROADCAST for fg processors is called in ERGASTULUM
664 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
665 time_Bcast=time_Bcast+MPI_Wtime()-time00
666 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
668 c call MPI_Barrier(FG_COMM,IERROR)
670 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
671 & MPI_DOUBLE_PRECISION,
672 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
677 c write (2,*) "My chunk of z"
678 c do i=1,3*my_ng_count
681 time_scatter=time_scatter+MPI_Wtime()-time00
683 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
691 temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
695 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
696 c write (2,*) "Before REDUCE"
697 c write (2,*) "d_a_tmp before reduce"
699 c write (2,*) i,temp(i)
703 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
704 & MPI_SUM,king,FG_COMM,IERR)
705 time_reduce=time_reduce+MPI_Wtime()-time00
706 c write (2,*) "After REDUCE"
718 d_a_tmp(ind)=d_a_tmp(ind)
719 & -fricmat(j,i)*z((j-1)*3+k+1)
724 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
729 c write (iout,*) "Vector d_a"
731 c write (2,*) i,d_a_tmp(i)