From: Adam Sieradzan Date: Mon, 18 Nov 2013 19:02:48 +0000 (+0100) Subject: dzialajacy gradient X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=39aaaf5f5ec4de449f18993aacdec4b43cabb1b0 dzialajacy gradient Conflicts: source/unres/src_MD-M/energy_p_new_barrier.F --- diff --git a/source/unres/src_MD-M/CMakeLists.txt b/source/unres/src_MD-M/CMakeLists.txt index cfdc654..3884bc1 100644 --- a/source/unres/src_MD-M/CMakeLists.txt +++ b/source/unres/src_MD-M/CMakeLists.txt @@ -203,7 +203,7 @@ set_property(SOURCE ${UNRES_MDM_SRC3} PROPERTY COMPILE_FLAGS ${FFLAGS3} ) #========================================= if(UNRES_MD_FF STREQUAL "GAB" ) # set preprocesor flags - set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" ) + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB -DNEWCORR" ) #========================================= # Settings for E0LL2Y force field diff --git a/source/unres/src_MD-M/COMMON.CONTACTS b/source/unres/src_MD-M/COMMON.CONTACTS index 4b81714..45c578b 100644 --- a/source/unres/src_MD-M/COMMON.CONTACTS +++ b/source/unres/src_MD-M/COMMON.CONTACTS @@ -28,7 +28,8 @@ C 10/30/99 Added other pre-computed vectors and matrices needed C to calculate three - six-order el-loc correlation terms double precision Ug,Ugder,Ug2,Ug2der,obrot,obrot2,obrot_der, & obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,gmu,gUb2 - & DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der + & DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der, + & gtEug common /rotat/ Ug(2,2,maxres),Ugder(2,2,maxres),Ug2(2,2,maxres), & Ug2der(2,2,maxres),obrot(2,maxres),obrot2(2,maxres), & obrot_der(2,maxres),obrot2_der(2,maxres) @@ -40,7 +41,7 @@ C amino-acid residue. & Dtobr2(2,maxres),Dtobr2der(2,maxres), & EUg(2,2,maxres),EUgder(2,2,maxres),CUg(2,2,maxres), & CUgder(2,2,maxres),DUg(2,2,maxres),Dugder(2,2,maxres), - & DtUg2(2,2,maxres),DtUg2der(2,2,maxres) + & DtUg2(2,2,maxres),DtUg2der(2,2,maxres),gtEUg(2,2,maxres) C This common block contains vectors and matrices dependent on two C consecutive amino-acid residues. double precision Ug2Db1t,Ug2Db1tder,CUgb2,CUgb2der,EUgC, diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index 95472be..3f1b981 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -25,16 +25,17 @@ C 6/23/01 - constants for double torsionals C 9/18/99 - added Fourier coeffficients of the expansion of local energy C surface double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde, - &bnew1,bnew2,eenew,ggb1,ggb2 + &bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee integer nloctyp common/fourier/ b1(2,maxres),b2(2,maxres), & bnew1(3,2,-maxtor:maxtor),bnew2(3,2,-maxtor:maxtor), & cc(2,2,-maxtor:maxtor), - & dd(2,2,-maxtor:maxtor),ee(2,2,-maxtor:maxtor), + & dd(2,2,-maxtor:maxtor),eeold(2,2,-maxtor:maxtor), & eenew(2,-maxtor:maxtor), + & ee(2,2,maxres), & ctilde(2,2,-maxtor:maxtor), & dtilde(2,2,-maxtor:maxtor),b1tilde(2,maxres), & b2tilde(2,maxres), - & gtb1(2,maxres),gtb2(2,maxres), + & gtb1(2,maxres),gtb2(2,maxres),gtEE(2,2,maxres), & nloctyp diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index c13d341..ab7d468 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2297,7 +2297,14 @@ c endif gtb1(2,i-2)=0.0 b2(2,i-2)=bnew2(1,2,iti) gtb2(2,i-2)=0.0 -c EE(1,1,iti)=0.0d0 + EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1)) + EE(1,2,i-2)=eeold(1,2,iti) + EE(2,1,i-2)=eeold(2,1,iti) + EE(2,2,i-2)=eeold(2,2,iti) + gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1)) + gtEE(1,2,i-2)=0.0d0 + gtEE(2,2,i-2)=0.0d0 + gtEE(2,1,i-2)=0.0d0 c EE(2,2,iti)=0.0d0 c EE(1,2,iti)=0.5d0*eenew(1,iti) c EE(2,1,iti)=0.5d0*eenew(1,iti) @@ -2383,7 +2390,6 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) Ug2der(2,2,i-2)=0.0d0 endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then -#ifndef NEWCORR if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) else @@ -2395,7 +2401,6 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then else iti1=ntortyp+1 endif -#endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) @@ -2407,7 +2412,13 @@ c if (i .gt. iatel_s+2) then call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" #endif - call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2)) +c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), +c & EE(1,2,iti),EE(2,2,iti) + call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) + call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) +c write(iout,*) "Macierz EUG", +c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2), +c & eug(2,2,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2)) @@ -2430,7 +2441,7 @@ c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" enddo endif call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2)) - call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2)) + call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) enddo @@ -2866,8 +2877,7 @@ C ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi num_conti=0 -c TU ZLE -c call eelecij(i,i+2,ees,evdw1,eel_loc) + call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) num_cont_hb(i)=num_conti enddo @@ -2885,8 +2895,8 @@ c call eelecij(i,i+2,ees,evdw1,eel_loc) ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi num_conti=num_cont_hb(i) -c TU ZLE -c call eelecij(i,i+3,ees,evdw1,eel_loc) +c write(iout,*) "JESTEM W PETLI" + call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti @@ -2894,8 +2904,8 @@ c call eelecij(i,i+3,ees,evdw1,eel_loc) c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c -c do i=iatel_s,iatel_e - do i=4,5 + do i=iatel_s,iatel_e +c do i=7,7 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2908,9 +2918,9 @@ c do i=iatel_s,iatel_e zmedi=c(3,i)+0.5d0*dzi c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) -c do j=ielstart(i),ielend(i) - do j=8,9 - write (iout,*) 'tu wchodze',i,j,itype(i),itype(j) + do j=ielstart(i),ielend(i) +c do j=13,13 +c write (iout,*) 'tu wchodze',i,j,itype(i),itype(j) if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -3640,7 +3650,9 @@ C Third- and fourth-order contributions from turns dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2), + & gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2), + & auxgmat2(2,2),auxgmatt2(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3664,9 +3676,24 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) +c auxalary matices for theta gradient +c auxalary matrix for i+1 and constant i+2 + call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1)) +c auxalary matrix for i+2 and constant i+1 + call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) + call transpose2(auxgmat1(1,1),auxgmatt1(1,1)) + call transpose2(auxgmat2(1,1),auxgmatt2(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) + call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) +C Derivatives in theta + gloc(nphi+i,icg)=gloc(nphi+i,icg) + & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 + gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) + & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) cd write (2,*) 'i,',i,' j',j,'eello_turn3', @@ -3740,7 +3767,11 @@ C Third- and fourth-order contributions from turns dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2), + & auxgEvec1(2),auxgEvec2(2),auxgEvec3(2), + & gte1t(2,2),gte2t(2,2),gte3t(2,2), + & gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2), + & gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3760,6 +3791,7 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn4(i,a_temp,eello_turn4_num) c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 +c write(iout,*)"WCHODZE W PROGRAM" a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -3771,37 +3803,83 @@ c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) call transpose2(Eug(1,1,i+3),e3t(1,1)) +C Ematrix derivative in theta + call transpose2(gtEUg(1,1,i+1),gte1t(1,1)) + call transpose2(gtEug(1,1,i+2),gte2t(1,1)) + call transpose2(gtEug(1,1,i+3),gte3t(1,1)) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) +c eta1 in derivative theta + call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) +c auxgvec is derivative of Ub2 so i+3 theta call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) +c auxalary matrix of E i+1 + call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1)) c s1=0.0 c gs1=0.0 s1=scalar2(b1(1,i+2),auxvec(1)) -c gs1=scalar2(gtb1(1,i+2),auxgvec(1)) +c derivative of theta i+2 with constant i+3 + gs23=scalar2(gtb1(1,i+2),auxvec(1)) +c derivative of theta i+2 with constant i+2 + gs32=scalar2(b1(1,i+2),auxgvec(1)) +c derivative of E matix in theta of i+1 + gsE13=scalar2(b1(1,i+2),auxgEvec1(1)) + call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) +c ea31 in derivative theta + call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) +c auxilary matrix auxgvec of Ub2 with constant E matirx call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) +c auxilary matrix auxgEvec1 of E matix with Ub2 constant + call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) + c s2=0.0 c gs2=0.0 s2=scalar2(b1(1,i+1),auxvec(1)) -c gs2=scalar2(gtb1(1,i+1),auxgvec(1)) -c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),ggb1(1,i+2), -c & ggb1(1,i+1) +c derivative of theta i+1 with constant i+3 + gs13=scalar2(gtb1(1,i+1),auxvec(1)) +c derivative of theta i+2 with constant i+1 + gs21=scalar2(b1(1,i+1),auxgvec(1)) +c derivative of theta i+3 with constant i+1 + gsE31=scalar2(b1(1,i+1),auxgEvec3(1)) +c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2), +c & gtb1(1,i+1) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) +c two derivatives over diffetent matrices +c gtae3e2 is derivative over i+3 + call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1)) +c ae3gte2 is derivative over i+2 + call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) +c three possible derivative over theta E matices +c i+1 + call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1)) +c i+2 + call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1)) +c i+3 + call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) + + gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2)) + gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2)) + gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2)) + eello_turn4=eello_turn4-(s1+s2+s3) #ifdef NEWCORR -c geel_loc_ij=-(gs1+gs2) -c gloc(nphi+i,icg)=gloc(nphi+i,icg)- -c & gs1 + gloc(nphi+i,icg)=gloc(nphi+i,icg) + & -(gs13+gsE13+gsEE1)*wturn4 + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg) + & -(gs23+gs21+gsEE2)*wturn4 + gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg) + & -(gs32+gsE31+gsEE3)*wturn4 c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- c & gs2 #endif if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn4',i,j,-(s1+s2+s3) -cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), -cd & ' eello_turn4_num',8*eello_turn4_num +c write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), +c & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) @@ -4544,7 +4622,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 - gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo C Ufff.... We've done all this!!! return @@ -4850,7 +4928,7 @@ c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) enddo return end diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 9db9b60..03efdb2 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -991,21 +991,22 @@ c Dtilde(1,1,i)=0.0d0 c Dtilde(1,2,i)=0.0d0 c Dtilde(2,1,i)=0.0d0 c Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= b(10)+b(11) - EE(2,2,i)=-b(10)+b(11) - EE(2,1,i)= b(12)-b(13) - EE(1,2,i)= b(12)+b(13) - EE(1,1,-i)= b(10)+b(11) - EE(2,2,-i)=-b(10)+b(11) - EE(2,1,-i)=-b(12)+b(13) - EE(1,2,-i)=-b(12)-b(13) - + EEold(1,1,i)= b(10)+b(11) + EEold(2,2,i)=-b(10)+b(11) + EEold(2,1,i)= b(12)-b(13) + EEold(1,2,i)= b(12)+b(13) + EEold(1,1,-i)= b(10)+b(11) + EEold(2,2,-i)=-b(10)+b(11) + EEold(2,1,-i)=-b(12)+b(13) + EEold(1,2,-i)=-b(12)-b(13) + write(iout,*) "TU DOCHODZE" c ee(1,1,i)=1.0d0 c ee(2,2,i)=1.0d0 c ee(2,1,i)=0.0d0 c ee(1,2,i)=0.0d0 c ee(2,1,i)=ee(1,2,i) enddo +c lprint=.true. if (lprint) then do i=1,nloctyp write (iout,*) 'Type',i @@ -1023,10 +1024,11 @@ c ee(2,1,i)=ee(1,2,i) enddo write(iout,*) 'EE' do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) + write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo endif +c lprint=.false. C C Read electrostatic-interaction parameters