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 .and. itype(i).ne.ntyp1) 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 .and. itype(i).ne.ntyp1) 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)
147 double precision work(8*maxres6)
148 integer iwork(maxres6)
149 common /przechowalnia/ Gcopy,Ghalf
151 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
152 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
154 c Determine the number of degrees of freedom (dimen) and the number of
156 dimen=(nct-nnt+1)+nside
157 dimen1=(nct-nnt)+(nct-nnt+1)
160 if (nfgtasks.gt.1) then
162 call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
163 time_Bcast=time_Bcast+MPI_Wtime()-time00
164 call int_bounds(dimen,igmult_start,igmult_end)
165 igmult_start=igmult_start-1
166 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
167 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
168 my_ng_count=igmult_end-igmult_start
169 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
170 & MPI_INTEGER,FG_COMM,IERROR)
171 write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
172 & ' absolute rank',myrank,' igmult_start',igmult_start,
173 & ' igmult_end',igmult_end,' count',my_ng_count
174 write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
175 write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
185 c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
186 c Zeroing out A and fricmat
192 c Diagonal elements of the dC part of A and the respective friction coefficients
204 c Off-diagonal elements of the dC part of A
211 c Diagonal elements of the dX part of A and the respective friction coefficients
221 massvec(ii)=msc(iabs(iti))
222 if (iti.ne.10 .and. iti.ne.ntyp1) then
226 Gmat(ii1,ii1)=ISC(iabs(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)
331 sqreig(i)=dsqrt(Geigen(i))
339 Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
340 Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
341 Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
346 write (iout,*) "Comparison of original and restored G"
349 write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
350 & Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
356 c-------------------------------------------------------------------------------
357 SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
358 implicit real*8 (a-h,o-z)
360 include 'COMMON.IOUNITS'
361 double precision A(LM2,LM3),B(LM2)
365 WRITE(IOUT,600) (I,I=KA,KB)
366 WRITE(IOUT,601) (B(I),I=KA,KB)
370 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
376 4 IF (KB.EQ.NC) RETURN
380 600 FORMAT (// 9H ROOT NO.,I4,9I11)
381 601 FORMAT (/5X,10(1PE11.4))
383 603 FORMAT (I5,10F11.5)
386 c-------------------------------------------------------------------------------
387 SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
388 implicit real*8 (a-h,o-z)
390 include 'COMMON.IOUNITS'
391 double precision A(LM2,LM3)
395 WRITE(IOUT,600) (I,I=KA,KB)
399 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
405 4 IF (KB.EQ.NC) RETURN
409 600 FORMAT (//5x,9I11)
411 603 FORMAT (I5,10F11.3)
414 c-------------------------------------------------------------------------------
415 SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
416 implicit real*8 (a-h,o-z)
418 include 'COMMON.IOUNITS'
419 double precision A(LM2,LM3)
423 WRITE(IOUT,600) (I,I=KA,KB)
427 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
433 4 IF (KB.EQ.NC) RETURN
437 600 FORMAT (//5x,7(3I5,2x))
439 603 FORMAT (I5,7(3F5.1,2x))
442 c-------------------------------------------------------------------------------
443 SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
444 implicit real*8 (a-h,o-z)
446 include 'COMMON.IOUNITS'
447 double precision A(LM2,LM3)
451 WRITE(IOUT,600) (I,I=KA,KB)
455 WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
461 4 IF (KB.EQ.NC) RETURN
465 600 FORMAT (//5x,4(3I9,2x))
467 603 FORMAT (I5,4(3F9.3,2x))
470 c---------------------------------------------------------------------------
471 SUBROUTINE ginv_mult(z,d_a_tmp)
472 implicit real*8 (a-h,o-z)
478 include 'COMMON.SETUP'
479 include 'COMMON.TIME1'
481 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
482 &,time01,zcopy(dimen3)
484 if (nfgtasks.gt.1) then
485 if (fg_rank.eq.0) then
486 c The matching BROADCAST for fg processors is called in ERGASTULUM
488 call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
489 time_Bcast=time_Bcast+MPI_Wtime()-time00
490 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
492 c write (2,*) "time00",time00
493 c write (2,*) "Before Scatterv"
495 c write (2,*) "Whole z (for FG master)"
499 c call MPI_Barrier(FG_COMM,IERROR)
501 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
502 & MPI_DOUBLE_PRECISION,
503 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
507 c write (2,*) "My chunk of z"
508 c do i=1,3*my_ng_count
511 c write (2,*) "After SCATTERV"
513 c write (2,*) "MPI_Wtime",MPI_Wtime()
514 time_scatter=time_scatter+MPI_Wtime()-time00
516 time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
518 c write (2,*) "time_scatter",time_scatter
519 c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
528 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
529 c & Ginv(i,j),z((j-1)*3+k+1),
530 c & Ginv(i,j)*z((j-1)*3+k+1)
531 c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
532 temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
536 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
537 c write (2,*) "Before REDUCE"
539 c write (2,*) "z before reduce"
541 c write (2,*) i,temp(i)
544 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
545 & MPI_SUM,king,FG_COMM,IERR)
546 time_reduce=time_reduce+MPI_Wtime()-time00
547 c write (2,*) "After REDUCE"
559 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
561 c & Ginv(i,j),z((j-1)*3+k+1),
562 c & Ginv(i,j)*z((j-1)*3+k+1)
563 d_a_tmp(ind)=d_a_tmp(ind)
564 & +Ginv(j,i)*z((j-1)*3+k+1)
565 c d_a_tmp(ind)=d_a_tmp(ind)
566 c & +Ginv(i,j)*z((j-1)*3+k+1)
571 time_ginvmult=time_ginvmult+MPI_Wtime()-time01
578 c---------------------------------------------------------------------------
580 SUBROUTINE ginv_mult_test(z,d_a_tmp)
583 c include 'COMMON.MD'
584 double precision z(dimen),d_a_tmp(dimen)
585 double precision ztmp(dimen/3),dtmp(dimen/3)
590 c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
599 ztmp(j)=z((j-1)*3+k+1)
602 call alignx(16,ztmp(1))
603 call alignx(16,dtmp(1))
604 call alignx(16,Ginv(1,1))
609 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
620 c---------------------------------------------------------------------------
621 SUBROUTINE fricmat_mult(z,d_a_tmp)
628 include 'COMMON.IOUNITS'
629 include 'COMMON.SETUP'
630 include 'COMMON.TIME1'
632 include 'COMMON.LANGEVIN'
634 include 'COMMON.LANGEVIN.lang0'
636 double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
637 &,time01,zcopy(dimen3)
639 if (nfgtasks.gt.1) then
640 if (fg_rank.eq.0) then
641 c The matching BROADCAST for fg processors is called in ERGASTULUM
643 call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
644 time_Bcast=time_Bcast+MPI_Wtime()-time00
645 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
647 c call MPI_Barrier(FG_COMM,IERROR)
649 call MPI_Scatterv(z,ng_counts(0),ng_start(0),
650 & MPI_DOUBLE_PRECISION,
651 & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
656 c write (2,*) "My chunk of z"
657 c do i=1,3*my_ng_count
660 time_scatter=time_scatter+MPI_Wtime()-time00
662 time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
670 temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
674 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
675 c write (2,*) "Before REDUCE"
676 c write (2,*) "d_a_tmp before reduce"
678 c write (2,*) i,temp(i)
682 call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
683 & MPI_SUM,king,FG_COMM,IERR)
684 time_reduce=time_reduce+MPI_Wtime()-time00
685 c write (2,*) "After REDUCE"
697 d_a_tmp(ind)=d_a_tmp(ind)
698 & -fricmat(j,i)*z((j-1)*3+k+1)
703 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
708 c write (iout,*) "Vector d_a"
710 c write (2,*) i,d_a_tmp(i)