From 2a20f5fb020ec764c5e3cd775cf7467fa75d22cc Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 14 Nov 2013 14:32:27 +0100 Subject: [PATCH] zabezpieczenie przedczarkowe --- PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new | 71 +++++ source/unres/src_MD-M/CMakeLists.txt | 2 +- source/unres/src_MD-M/COMMON.CONTACTS | 3 +- source/unres/src_MD-M/COMMON.TORSION | 15 +- source/unres/src_MD-M/energy_p_new_barrier.F | 357 ++++++++++++++++-------- source/unres/src_MD-M/parmread.F | 32 ++- 6 files changed, 349 insertions(+), 131 deletions(-) create mode 100644 PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new diff --git a/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new b/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new new file mode 100644 index 0000000..231618e --- /dev/null +++ b/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new @@ -0,0 +1,71 @@ + 3 # Number of local interaction types +Gly + 0.000000000000000 + 0.791965124028570 + 0.206068961118571 + 0.000000000000000 + 0.000000000000000 + 2.373462483972307 + -0.927962753087420 + 0.000000000000000 + 0.000000000000000 + 1.329421814829764 + -0.370576187607876 + 0.000000000000000 + 0.000000000000000 + 1.000000000000000 + 0.500000000000000 + 0.500000000000000 + 1.000000000000000 + 1.000000000000000 + 1.500000000000000 + 1.500000000000000 + 2.000000000000000 + 2.000000000000000 +Ala + 0.000000000000000 + 0.500261572719827 + -0.233786079150650 + -0.878020534542259 + 1.501220349138902 + -2.089734050079038 + 2.302365702770721 + -0.532502145482045 + -1.596421505165690 + 1.276301651241011 + 0.399942874780603 + -0.543778421330412 + 0.400478498916355 + 1.000000000000000 + 0.500000000000000 + 0.500000000000000 + 1.000000000000000 + 1.000000000000000 + 1.500000000000000 + 1.500000000000000 + 2.000000000000000 + 2.000000000000000 +Pro + 0.000000000000000 + -1.286480000000000 + 0.031808800000000 + -0.906628000000000 + 1.015390000000000 + -1.548870000000000 + 1.914370000000000 + 0.664143000000000 + -0.454839000000000 + -0.051291400000000 + 0.103179000000000 + 0.316367000000000 + 0.045277200000000 + 1.000000000000000 + 0.500000000000000 + 0.500000000000000 + 1.000000000000000 + 1.000000000000000 + 1.500000000000000 + 1.500000000000000 + 2.000000000000000 + 2.000000000000000 + diff --git a/source/unres/src_MD-M/CMakeLists.txt b/source/unres/src_MD-M/CMakeLists.txt index cba4cc2..cfdc654 100644 --- a/source/unres/src_MD-M/CMakeLists.txt +++ b/source/unres/src_MD-M/CMakeLists.txt @@ -210,7 +210,7 @@ if(UNRES_MD_FF STREQUAL "GAB" ) #========================================= elseif(UNRES_MD_FF STREQUAL "E0LL2Y") # set preprocesor flags - set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB" ) + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB -DNEWCORR" ) endif(UNRES_MD_FF STREQUAL "GAB") diff --git a/source/unres/src_MD-M/COMMON.CONTACTS b/source/unres/src_MD-M/COMMON.CONTACTS index 5b3a90d..4b81714 100644 --- a/source/unres/src_MD-M/COMMON.CONTACTS +++ b/source/unres/src_MD-M/COMMON.CONTACTS @@ -27,7 +27,7 @@ c & dipderx(3,5,4,maxconts,maxres) 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, + & obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,gmu,gUb2 & DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der common /rotat/ Ug(2,2,maxres),Ugder(2,2,maxres),Ug2(2,2,maxres), & Ug2der(2,2,maxres),obrot(2,maxres),obrot2(2,maxres), @@ -35,6 +35,7 @@ C to calculate three - six-order el-loc correlation terms C This common block contains vectors and matrices dependent on a single C amino-acid residue. common /precomp1/ mu(2,maxres),muder(2,maxres),Ub2(2,maxres), + & gmu(2,maxres),gUb2(2,maxres), & Ub2der(2,maxres),Ctobr(2,maxres),Ctobrder(2,maxres), & Dtobr2(2,maxres),Dtobr2der(2,maxres), & EUg(2,2,maxres),EUgder(2,2,maxres),CUg(2,2,maxres), diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index 5f4b2b7..95472be 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -24,10 +24,17 @@ C 6/23/01 - constants for double torsionals & ntermd_2(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2) 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 + double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde, + &bnew1,bnew2,eenew,ggb1,ggb2 integer nloctyp - common/fourier/ b1(2,-maxtor:maxtor),b2(2,-maxtor:maxtor) - & ,cc(2,2,-maxtor:maxtor), + 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), + & eenew(2,-maxtor:maxtor), & ctilde(2,2,-maxtor:maxtor), - & dtilde(2,2,-maxtor:maxtor),b1tilde(2,-maxtor:maxtor),nloctyp + & dtilde(2,2,-maxtor:maxtor),b1tilde(2,maxres), + & b2tilde(2,maxres), + & gtb1(2,maxres),gtb2(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 d4e2092..c13d341 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2252,11 +2252,70 @@ C C Compute the virtual-bond-torsional-angle dependent quantities needed C to calculate the el-loc multibody terms of various order. C + write(iout,*) 'nphi=',nphi,nres #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else do i=3,nres+1 #endif +#ifdef NEWCORR + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itortyp(itype(i-2)) + else + iti=ntortyp+1 + endif +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itortyp(itype(i-1)) + else + iti1=ntortyp+1 + endif + write(iout,*),i + b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0) + & +bnew1(2,1,iti)*sin(theta(i-1)) + & +bnew1(3,1,iti)*cos(theta(i-1)/2.0) + gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0 + & +bnew1(2,1,iti)*cos(theta(i-1)) + & -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0 +c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) +c &*(cos(theta(i)/2.0) + b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0) + & +bnew2(2,1,iti)*sin(theta(i-1)) + & +bnew2(3,1,iti)*cos(theta(i-1)/2.0) +c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) +c &*(cos(theta(i)/2.0) + gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0 + & +bnew2(2,1,iti)*cos(theta(i-1)) + & -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0 +c if (ggb1(1,i).eq.0.0d0) then +c write(iout,*) 'i=',i,ggb1(1,i), +c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0, +c &bnew1(2,1,iti)*cos(theta(i)), +c &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0 +c endif + b1(2,i-2)=bnew1(1,2,iti) + 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 +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) +c b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i)) +c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) + b1tilde(1,i-2)=b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)=b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) + write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) + write (iout,*) 'theta=', theta(i-1) + enddo +#ifdef PARMAT + do i=ivec_start+2,ivec_end+2 +#else + do i=3,nres+1 +#endif +#endif if (i .lt. nres+1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) @@ -2324,6 +2383,7 @@ C 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 @@ -2335,13 +2395,18 @@ 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) cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then - call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2)) + call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) +#ifdef NEWCORR + 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)) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then @@ -2364,7 +2429,7 @@ c if (i .gt. iatel_s+2) then enddo enddo endif - call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2)) + 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)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) @@ -2380,7 +2445,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then iti1=ntortyp+1 endif do k=1,2 - mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) + mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo cd write (iout,*) 'mu ',mu(:,i-2) cd write (iout,*) 'mu1',mu1(:,i-2) @@ -2393,7 +2458,7 @@ cd write (iout,*) 'mu2',mu2(:,i-2) call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) C Vectors and matrices dependent on a single virtual-bond dihedral. - call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) + call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2)) @@ -2709,7 +2774,7 @@ C dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -2801,7 +2866,8 @@ C ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi num_conti=0 - call eelecij(i,i+2,ees,evdw1,eel_loc) +c TU ZLE +c 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 @@ -2819,7 +2885,8 @@ C ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi num_conti=num_cont_hb(i) - call eelecij(i,i+3,ees,evdw1,eel_loc) +c TU ZLE +c 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 @@ -2827,7 +2894,8 @@ C c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c - do i=iatel_s,iatel_e +c do i=iatel_s,iatel_e + do i=4,5 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2840,8 +2908,9 @@ c zmedi=c(3,i)+0.5d0*dzi c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) - do j=ielstart(i),ielend(i) -c write (iout,*) i,j,itype(i),itype(j) +c do j=ielstart(i),ielend(i) + do j=8,9 + 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 @@ -2880,7 +2949,8 @@ C------------------------------------------------------------------------------- dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), + & gmuij2(4),gmuji2(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -3097,6 +3167,7 @@ C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al. C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms C are computed for EVERY pair of non-contiguous peptide groups. C + if (j.lt.nres-1) then j1=j+1 j2=j-1 @@ -3105,10 +3176,19 @@ C j2=j-2 endif kkk=0 + lll=0 do k=1,2 do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +#ifdef NEWCORR + gmuij1(kkk)=gtb1(k,i)*mu(l,j) + write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i + gmuij2(kkk)=gUb2(k,i-1)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j) + write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j + gmuji2(kkk)=mu(k,i)*gUb2(l,j-1) +#endif enddo enddo cd write (iout,*) 'EELEC: i',i,' j',j @@ -3274,6 +3354,40 @@ cgrad endif C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) +C Calculate patrial derivative for theta angle +#ifdef NEWCORR + geel_loc_ij=a22*gmuij1(1) + & +a23*gmuij1(2) + & +a32*gmuij1(3) + & +a33*gmuij1(4) + write(iout,*) "derivative over thatai" + write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), + & a33*gmuij1(4) + gloc(nphi+i,icg)=gloc(nphi+i,icg)+ + & geel_loc_ij*wel_loc + write(iout,*) "derivative over thatai-1" + write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), + & a33*gmuij2(4) + geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3) + & +a33*gmuij2(4) + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ + & geel_loc_ij*wel_loc + geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3) + & +a33*gmuji1(4) + write(iout,*) "derivative over thataj" + write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), + & a33*gmuji1(4) + + gloc(nphi+j,icg)=gloc(nphi+j,icg)+ + & geel_loc_ji*wel_loc + geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3) + & +a33*gmuji2(4) + write(iout,*) "derivative over thataj-1" + write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), + & a33*gmuji2(4) + gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ + & geel_loc_ji*wel_loc +#endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') @@ -3626,7 +3740,7 @@ 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),auxgvec(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, @@ -3659,14 +3773,31 @@ c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(Eug(1,1,i+3),e3t(1,1)) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(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)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(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) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(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 +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), @@ -3675,14 +3806,14 @@ 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)) call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1)) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3690,10 +3821,10 @@ C Derivatives in gamma(i+1) C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1)) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3708,10 +3839,10 @@ C Derivatives of this turn contributions in DC(i+2) a_temp(2,2)=agg(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3727,10 +3858,10 @@ C Remaining derivatives of this turn contribution a_temp(2,2)=aggi(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3741,10 +3872,10 @@ C Remaining derivatives of this turn contribution a_temp(2,2)=aggi1(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3755,10 +3886,10 @@ C Remaining derivatives of this turn contribution a_temp(2,2)=aggj(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -3769,10 +3900,10 @@ C Remaining derivatives of this turn contribution a_temp(2,2)=aggj1(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) @@ -6816,10 +6947,10 @@ C--------------------------------------------------------------------------- do iii=1,2 dipi(iii,1)=Ub2(iii,i) dipderi(iii)=Ub2der(iii,i) - dipi(iii,2)=b1(iii,iti1) + dipi(iii,2)=b1(iii,i+1) dipj(iii,1)=Ub2(iii,j) dipderj(iii)=Ub2der(iii,j) - dipj(iii,2)=b1(iii,itj1) + dipj(iii,2)=b1(iii,j+1) enddo kkk=0 do iii=1,2 @@ -6999,26 +7130,26 @@ C They are needed only when the fifth- or the sixth-order cumulants are C indluded. IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN call transpose2(AEA(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) call transpose2(AEAderg(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) - call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) - call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1)) call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) call transpose2(AEA(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2)) call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2)) call transpose2(AEAderg(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2)) - call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2)) - call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2)) call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2)) call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2)) call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2)) @@ -7027,20 +7158,20 @@ C Calculate the Cartesian derivatives of the vectors. do kkk=1,5 do lll=1,3 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti), + call matvec2(auxmat(1,1),b1(1,i), & AEAb1derx(1,lll,kkk,iii,1,1)) call matvec2(auxmat(1,1),Ub2(1,i), & AEAb2derx(1,lll,kkk,iii,1,1)) - call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & AEAb1derx(1,lll,kkk,iii,2,1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), & AEAb2derx(1,lll,kkk,iii,2,1)) call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj), + call matvec2(auxmat(1,1),b1(1,j), & AEAb1derx(1,lll,kkk,iii,1,2)) call matvec2(auxmat(1,1),Ub2(1,j), & AEAb2derx(1,lll,kkk,iii,1,2)) - call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1), & AEAb1derx(1,lll,kkk,iii,2,2)) call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1), & AEAb2derx(1,lll,kkk,iii,2,2)) @@ -7137,26 +7268,26 @@ C indluded. IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or. & (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN call transpose2(AEA(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) call transpose2(AEAderg(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) - call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) - call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1)) call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) call transpose2(AEA(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2)) call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2)) call transpose2(AEAderg(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2)) - call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2)) - call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2)) call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2)) call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2)) call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2)) @@ -7165,20 +7296,20 @@ C Calculate the Cartesian derivatives of the vectors. do kkk=1,5 do lll=1,3 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti), + call matvec2(auxmat(1,1),b1(1,i), & AEAb1derx(1,lll,kkk,iii,1,1)) call matvec2(auxmat(1,1),Ub2(1,i), & AEAb2derx(1,lll,kkk,iii,1,1)) - call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & AEAb1derx(1,lll,kkk,iii,2,1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), & AEAb2derx(1,lll,kkk,iii,2,1)) call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itl), + call matvec2(auxmat(1,1),b1(1,l), & AEAb1derx(1,lll,kkk,iii,1,2)) call matvec2(auxmat(1,1),Ub2(1,l), & AEAb2derx(1,lll,kkk,iii,1,2)) - call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1), + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1), & AEAb1derx(1,lll,kkk,iii,2,2)) call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j), & AEAb2derx(1,lll,kkk,iii,2,2)) @@ -7475,7 +7606,7 @@ C Contribution from graph II call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) - eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk)) + eello5_2=scalar2(AEAb1(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(k-1)=g_corr5_loc(k-1) @@ -7485,11 +7616,11 @@ C Explicit gradient in virtual-dihedral angles. vv(2)=pizda(2,1)-pizda(1,2) if (l.eq.j+1) then g_corr5_loc(l-1)=g_corr5_loc(l-1) - & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k))) else g_corr5_loc(j-1)=g_corr5_loc(j-1) - & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k))) endif C Cartesian gradient @@ -7501,7 +7632,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,iii)=derx(lll,kkk,iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k)) enddo enddo @@ -7556,7 +7687,7 @@ cd1110 continue call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) - eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl)) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(l-1)=g_corr5_loc(l-1) @@ -7565,7 +7696,7 @@ C Explicit gradient in virtual-dihedral angles. vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) g_corr5_loc(k-1)=g_corr5_loc(k-1) - & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl)) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l))) C Cartesian gradient do iii=1,2 @@ -7576,7 +7707,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,iii)=derx(lll,kkk,iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l)) enddo enddo @@ -7629,7 +7760,7 @@ C Contribution from graph IV call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) - eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj)) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(j-1)=g_corr5_loc(j-1) @@ -7638,7 +7769,7 @@ C Explicit gradient in virtual-dihedral angles. vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) g_corr5_loc(k-1)=g_corr5_loc(k-1) - & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj)) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j))) C Cartesian gradient do iii=1,2 @@ -7649,7 +7780,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j)) enddo enddo @@ -7931,8 +8062,8 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) - vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk) - vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk) + vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k) + vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k) s5=scalar2(vv(1),Dtobr2(1,i)) cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5) @@ -7945,8 +8076,8 @@ cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1)) vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) - vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk) - vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk) + vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k) + vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k) if (l.eq.j+1) then g_corr6_loc(l-1)=g_corr6_loc(l-1) & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i)) @@ -7985,10 +8116,10 @@ cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) - vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk) - & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk) - vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk) - & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk) + vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k) + & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k) + vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k) + & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k) s5=scalar2(vv(1),Dtobr2(1,i)) derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5) enddo @@ -8228,10 +8359,10 @@ C energy moment and not to the cluster cumulant. #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) #endif - call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) - call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1)) + s2=0.5d0*scalar2(b1(1,k),auxvec(1)) + call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) call transpose2(EE(1,1,itk),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) @@ -8246,13 +8377,13 @@ cd & "sum",-(s2+s3+s4) #endif c eello6_graph3=-s4 C Derivatives in gamma(k-1) - call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k)) g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4) C Derivatives in gamma(l-1) - call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) + call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1)) + s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8269,12 +8400,12 @@ C Cartesian derivatives. s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k) endif #endif - call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) - call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + s2=0.5d0*scalar2(b1(1,k),auxvec(1)) + call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1), & auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1), & pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) @@ -8362,11 +8493,11 @@ cd & ' itl',itl,' itl1',itl1 call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1)) s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) if (j.eq.l+1) then - call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif call transpose2(EUg(1,1,k),auxmat(1,1)) call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1)) @@ -8390,11 +8521,11 @@ C Derivatives in gamma(i-1) #endif s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1)) if (j.eq.l+1) then - call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i)) if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then @@ -8423,11 +8554,11 @@ C Derivatives in gamma(k-1) call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1)) s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1)) if (j.eq.l+1) then - call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif call transpose2(EUgder(1,1,k),auxmat1(1,1)) call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1)) @@ -8493,12 +8624,12 @@ C Cartesian derivatives. s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) if (j.eq.l+1) then call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), - & b1(1,itj1),auxvec(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec(1)) + & b1(1,j+1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec(1)) else call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), - & b1(1,itl1),auxvec(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec(1)) + & b1(1,l+1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec(1)) endif call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1), & pizda(1,1)) @@ -8598,12 +8729,12 @@ cd write (2,*) 'eello6_5',eello6_5 #ifdef MOMENT call transpose2(AEA(1,1,1),auxmat(1,1)) call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1)) - ss1=scalar2(Ub2(1,i+2),b1(1,itl)) + ss1=scalar2(Ub2(1,i+2),b1(1,l)) s1 = (auxmat(1,1)+auxmat(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1)) call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1)) - s2 = scalar2(b1(1,itk),vtemp1(1)) + s2 = scalar2(b1(1,k),vtemp1(1)) #ifdef MOMENT call transpose2(AEA(1,1,2),atemp(1,1)) call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1)) @@ -8618,7 +8749,7 @@ cd write (2,*) 'eello6_5',eello6_5 call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1)) call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1)) call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1)) - ss13 = scalar2(b1(1,itk),vtemp4(1)) + ss13 = scalar2(b1(1,k),vtemp4(1)) s13 = (gtemp(1,1)+gtemp(2,2))*ss13 #endif c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13 @@ -8652,12 +8783,12 @@ C Derivatives in gamma(i+3) #ifdef MOMENT call transpose2(AEA(1,1,1),auxmatd(1,1)) call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) - ss1d=scalar2(Ub2der(1,i+2),b1(1,itl)) + ss1d=scalar2(Ub2der(1,i+2),b1(1,l)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d #endif - call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1)) call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1)) s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1)) @@ -8705,9 +8836,9 @@ C Derivatives in gamma(i+5) call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1)) call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call transpose2(AEA(1,1,2),atempd(1,1)) call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1)) @@ -8717,7 +8848,7 @@ C Derivatives in gamma(i+5) s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) #ifdef MOMENT call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1)) - ss13d = scalar2(b1(1,itk),vtemp4d(1)) + ss13d = scalar2(b1(1,k),vtemp4d(1)) s13d = (gtemp(1,1)+gtemp(2,2))*ss13d #endif c s1d=0.0d0 @@ -8741,10 +8872,10 @@ C Cartesian derivatives call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1), & vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1)) call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1)) @@ -8788,7 +8919,7 @@ c s13d=0.0d0 derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4), & vtemp4d(1)) - ss13d = scalar2(b1(1,itk),vtemp4d(1)) + ss13d = scalar2(b1(1,k),vtemp4d(1)) s13d = (gtemp(1,1)+gtemp(2,2))*ss13d derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d enddo diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 6ea0840..9db9b60 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -907,29 +907,37 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp + do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) +#ifdef NEWCORR + read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1) +#endif if (lprint) then write (iout,*) 'Type',i write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13) endif - B1(1,i) = b(3) - B1(2,i) = b(5) - B1(1,-i) = b(3) - B1(2,-i) = -b(5) +c B1(1,i) = b(3) +c B1(2,i) = b(5) +c B1(1,-i) = b(3) +c B1(2,-i) = -b(5) c b1(1,i)=0.0d0 c b1(2,i)=0.0d0 - B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) +c B1tilde(1,i) = b(3) +c B1tilde(2,i) =-b(5) +c B1tilde(1,-i) =-b(3) +c B1tilde(2,-i) =b(5) c b1tilde(1,i)=0.0d0 c b1tilde(2,i)=0.0d0 - B2(1,i) = b(2) - B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) +c B2(1,i) = b(2) +c B2(2,i) = b(4) +c B2(1,-i) =b(2) +c B2(2,-i) =-b(4) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 -- 1.7.9.5