X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=21ac2a76528f88ff241c79c7b0d967aa1987b09d;hb=c50460147146c32657710433c7464783b5347bf5;hp=51fbd11688a474e0816b150a8b700403533ed3c2;hpb=478a9d9a1c99eb3f4bc4ca676ff3162bdd01d633;p=unres.git diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 51fbd11..21ac2a7 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -471,7 +471,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'mpif.h' #endif double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres) + & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' @@ -483,6 +483,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.TIME1' include 'COMMON.MAXGRAD' + include 'COMMON.SCCOR' #ifdef TIMING #ifdef MPI time01=MPI_Wtime() @@ -493,8 +494,9 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef DEBUG write (iout,*) "sum_gradient gvdwc, gvdwx" do i=1,nres - write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') - & i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3) + write (iout,'(i3,3f10.5,5x,3f10.5,5x,3f10.5,5x,3f10.5)') + & i,(gvdwx(j,i),j=1,3),(gvdwcT(j,i),j=1,3),(gvdwc(j,i),j=1,3), + & (gvdwcT(j,i),j=1,3) enddo call flush(iout) #endif @@ -527,6 +529,21 @@ c enddo call flush(iout) #endif #ifdef SPLITELE +#ifdef TSCSC + do i=1,nct + do j=1,3 + gradbufc(j,i)=wsc*gvdwc(j,i)+wsc*wscT*gvdwcT(j,i)+ + & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ + & welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+ + & wel_loc*gel_loc_long(j,i)+ + & wcorr*gradcorr_long(j,i)+ + & wcorr5*gradcorr5_long(j,i)+ + & wcorr6*gradcorr6_long(j,i)+ + & wturn6*gcorr6_turn_long(j,i)+ + & wstrain*ghpbc(j,i) + enddo + enddo +#else do i=1,nct do j=1,3 gradbufc(j,i)=wsc*gvdwc(j,i)+ @@ -540,6 +557,7 @@ c enddo & wstrain*ghpbc(j,i) enddo enddo +#endif #else do i=1,nct do j=1,3 @@ -571,6 +589,16 @@ c enddo gradbufc_sum(j,i)=gradbufc(j,i) enddo enddo +c call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres, +c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR) +c time_reduce=time_reduce+MPI_Wtime()-time00 +#ifdef DEBUG +c write (iout,*) "gradbufc_sum after allreduce" +c do i=1,nres +c write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3) +c enddo +c call flush(iout) +#endif #ifdef TIMING c time_allreduce=time_allreduce+MPI_Wtime()-time00 #endif @@ -585,6 +613,17 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 & " jgrad_end ",jgrad_end(i), & i=igrad_start,igrad_end) #endif +c +c Obsolete and inefficient code; we can make the effort O(n) and, therefore, +c do not parallelize this part. +c +c do i=igrad_start,igrad_end +c do j=jgrad_start(i),jgrad_end(i) +c do k=1,3 +c gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j) +c enddo +c enddo +c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo @@ -623,6 +662,16 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo enddo +c do i=nnt,nres-1 +c do k=1,3 +c gradbufc(k,i)=0.0d0 +c enddo +c do j=i+1,nres +c do k=1,3 +c gradbufc(k,i)=gradbufc(k,i)+gradbufc(k,j) +c enddo +c enddo +c enddo #ifdef DEBUG write (iout,*) "gradbufc after summing" do i=1,nres @@ -661,7 +710,7 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ - & welec*gelc_long(j,i) + & welec*gelc_long(j,i)+ & wel_loc*gel_loc_long(j,i)+ & wcorr*gcorr_long(j,i)+ & wcorr5*gradcorr5_long(j,i)+ @@ -677,11 +726,20 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) #endif +#ifdef TSCSC + gradx(j,i,icg)=wsc*gvdwx(j,i)+wsc*wscT*gvdwxT(j,i)+ + & wscp*gradx_scp(j,i)+ + & wbond*gradbx(j,i)+ + & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ + & wsccor*gsccorx(j,i) + & +wscloc*gsclocx(j,i) +#else gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) +#endif enddo enddo #ifdef DEBUG @@ -698,7 +756,6 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 & +wturn3*gel_loc_turn3(i) & +wturn6*gel_loc_turn6(i) & +wel_loc*gel_loc_loc(i) - & +wsccor*gsccor_loc(i) enddo #ifdef DEBUG write (iout,*) "gloc after adding corr" @@ -717,6 +774,19 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo +#ifdef DEBUG + write (iout,*) "gloc_sc before reduce" + do i=1,nres + do j=1,3 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif + do i=1,nres + do j=1,3 + gloc_scbuf(j,i)=gloc_sc(j,i,icg) + enddo + enddo time00=MPI_Wtime() call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 @@ -727,8 +797,18 @@ c time_allreduce=time_allreduce+MPI_Wtime()-time00 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 #ifdef DEBUG + write (iout,*) "gloc_sc after reduce" + do i=1,nres + do j=1,3 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif +#ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres write (iout,*) i,gloc(i,icg) @@ -764,6 +844,10 @@ c do i=1,nct gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i))) if (gvdwc_norm.gt.gvdwc_max) gvdwc_max=gvdwc_norm +#ifdef TSCSC + gvdwc_norm=dsqrt(scalar(gvdwcT(1,i),gvdwcT(1,i))) + if (gvdwc_norm.gt.gvdwc_max) gvdwc_max=gvdwc_norm +#endif gvdwc_scp_norm=dsqrt(scalar(gvdwc_scp(1,i),gvdwc_scp(1,i))) if (gvdwc_scp_norm.gt.gvdwc_scp_max) & gvdwc_scp_max=gvdwc_scp_norm @@ -802,6 +886,10 @@ c if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i))) if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm +#ifdef TSCSC + gvdwx_norm=dsqrt(scalar(gvdwxT(1,i),gvdwxT(1,i))) + if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm +#endif gradx_scp_norm=dsqrt(scalar(gradx_scp(1,i),gradx_scp(1,i))) if (gradx_scp_norm.gt.gradx_scp_max) & gradx_scp_max=gradx_scp_norm @@ -973,7 +1061,7 @@ C------------------------------------------------------------------------ & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ - & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, + & 'EHPB= ',1pE16.6,' WEIGHT=',1pD16.6, & ' (SS bridges & dist. cnstr.)'/ & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ @@ -4174,49 +4262,90 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij cd write (iout,*) "eij",eij + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - fac=waga*rdis/dd + fac=waga*rdis/dd + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo endif enddo ehpb=0.5D0*ehpb @@ -5585,7 +5714,7 @@ C Proline-Proline pair is a special case... & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci -c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) + write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo ! 6/20/98 - dihedral angle constraints edihcnstr=0.0d0 @@ -5700,6 +5829,7 @@ c do i=1,ndih_constr else difi=0.0 endif +c write (iout,*) "gloci", gloc(i-3,icg) cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, cd & rad2deg*phi0(i), rad2deg*drange(i), cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -5736,7 +5866,6 @@ c lprn=.true. phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 -C Regular cosine and sine terms do j=1,ntermd_1(itori,itori1,itori2) v1cij=v1c(1,j,itori,itori1,itori2) v1sij=v1s(1,j,itori,itori1,itori2) @@ -5771,6 +5900,7 @@ C Regular cosine and sine terms enddo gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 +c write (iout,*) "gloci", gloc(i-3,icg) enddo return end @@ -5803,26 +5933,57 @@ C Set lprn=.true. for debugging c lprn=.true. c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 - do i=iphi_start,iphi_end + do i=itau_start,itau_end esccor_ii=0.0D0 - itori=itype(i-2) - itori1=itype(i-1) + isccori=isccortyp(itype(i-2)) + isccori1=isccortyp(itype(i-1)) phii=phi(i) +cccc Added 9 May 2012 +cc Tauangle is torsional engle depending on the value of first digit +c(see comment below) +cc Omicron is flat angle depending on the value of first digit +c(see comment below) + + + do intertyp=1,3 !intertyp +cc Added 09 May 2012 (Adasko) +cc Intertyp means interaction type of backbone mainchain correlation: +c 1 = SC...Ca...Ca...Ca +c 2 = Ca...Ca...Ca...SC +c 3 = SC...Ca...Ca...SCi gloci=0.0D0 - do j=1,nterm_sccor - v1ij=v1sccor(j,itori,itori1) - v2ij=v2sccor(j,itori,itori1) - cosphi=dcos(j*phii) - sinphi=dsin(j*phii) + if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. + & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. + & (itype(i-1).eq.21))) + & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) + & .or.(itype(i-2).eq.21))) + & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. + & (itype(i-1).eq.21)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) + & cycle + do j=1,nterm_sccor(isccori,isccori1) + v1ij=v1sccor(j,intertyp,isccori,isccori1) + v2ij=v2sccor(j,intertyp,isccori,isccori1) + cosphi=dcos(j*tauangle(intertyp,i)) + sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo + gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci +c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi +c &gloc_sc(intertyp,i-3,icg) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6) + & (v1sccor(j,intertyp,itori,itori1),j=1,6) + & ,(v2sccor(j,intertyp,itori,itori1),j=1,6) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci + enddo !intertyp enddo +c do i=1,nres +c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg) +c enddo return end c---------------------------------------------------------------------------- @@ -7929,9 +8090,9 @@ C C Parallel Antiparallel C C o o -C /l\ /j\ -C / \ / \ -C /| o | | o |\ +C /l\ /j\ +C / \ / \ +C /| o | | o |\ C \ j|/k\| / \ |/k\|l / C \ / \ / \ / \ / C o o o o @@ -8026,22 +8187,22 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C \ /l\ /j\ / -C \ / \ / \ / -C o| o | | o |o -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C \ /l\ /j\ / C +C \ / \ / \ / C +C o| o | | o |o C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8212,18 +8373,18 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C j|/k\| / |/k\|l / -C / \ / / \ / -C / o / o -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C j|/k\| / |/k\|l / C +C / \ / / \ / C +C / o / o C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective @@ -8329,18 +8490,18 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C -C Parallel Antiparallel -C -C o o -C /l\ / \ /j\ -C / \ / \ / \ -C /| o |o o| o |\ -C \ j|/k\| \ |/k\|l -C \ / \ \ / \ -C o \ o \ -C i i -C +C C +C Parallel Antiparallel C +C C +C o o C +C /l\ / \ /j\ C +C / \ / \ / \ C +C /| o |o o| o |\ C +C \ j|/k\| \ |/k\|l C +C \ / \ \ / \ C +C o \ o \ C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective