X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=730500e14a409aed49325fa215f0c2bb0f3693f8;hb=4ed0b0657cdb9599bd714a010134061c94f509cc;hp=eb6146af1b6eb7f56d339dbd8535276c512d73de;hpb=05c18ccfcefc037cb8394bd779a4e51b7cbea6ec;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index eb6146a..730500e 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -224,6 +224,11 @@ c detecting NaNQ #ifdef MPL c endif #endif +C#define DEBUG +#ifdef DEBUG + call enerprint(energia,fact) +#endif +C#undef DEBUG if (calc_grad) then C C Sum up the components of the Cartesian gradient. @@ -265,12 +270,29 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i) & +fact(1)*wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + endif enddo @@ -309,12 +331,29 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+ & fact(1)*wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + endif enddo #endif @@ -895,7 +934,7 @@ C include 'COMMON.SBRIDGE' logical lprn common /srutu/icall - integer icant + integer icant,xshift,yshift,zshift external icant do i=1,210 do j=1,2 @@ -1034,6 +1073,13 @@ C lipbufthick is thickenes of lipid buffore & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 +C if (aa.ne.aa_aq(itypi,itypj)) then + +C write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa, +C & bb_aq(itypi,itypj)-bb, +C & sslipi,sslipj +C endif + C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj) C checking the distance dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 @@ -1116,6 +1162,7 @@ c & aux*e2/eps(itypi,itypj) c if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa +C#define DEBUG #ifdef DEBUG write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, @@ -1125,6 +1172,7 @@ c if (lprn) then & evdwij write (iout,*) "partial sum", evdw, evdw_t #endif +C#undef DEBUG c endif if (calc_grad) then C Calculate gradient components. @@ -2098,15 +2146,15 @@ C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (i.eq.1) then +C if (i.eq.1) then if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1) cycle - else - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1) cycle +C else +C if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 &) cycle - endif +C endif if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -2126,16 +2174,16 @@ C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e num_conti=0 C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (j.eq.1) then - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - &) cycle - else + if (j.le.1) cycle +C if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C &) cycle +C else if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C & .or.itype(j-1).eq.ntyp1 &) cycle - endif +C endif C C) cycle if (itel(j).eq.0) goto 1216 @@ -2225,6 +2273,12 @@ c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) if (shield_mode.gt.0) then +C#define DEBUG +#ifdef DEBUG + write(iout,*) "ees_compon",i,j,el1,el2, + & fac_shield(i),fac_shield(j) +#endif +C#undef DEBUG C fac_shield(i)=0.4 C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 @@ -2960,6 +3014,8 @@ C Derivatives due to the contact function & *fac_shield(i)*fac_shield(j) gacontp_hb3(k,num_conti,i)=gggp(k) + & *fac_shield(i)*fac_shield(j) + gacontm_hb1(k,num_conti,i)=ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) @@ -3028,6 +3084,18 @@ C Third- and fourth-order contributions from turns & aggj(3,4),aggj1(3,4),a_temp(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 if (j.eq.i+2) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) +C end of changes suggested by Ana + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 179 + CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Third-order contributions @@ -3107,6 +3175,7 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),pizda(1,1)) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) @@ -3155,17 +3224,19 @@ C Cartesian derivatives enddo endif + 179 continue else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+5).gt.nres) - & .or.((i-1).le.0) +C & .or.((i+5).gt.nres) +C & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1) goto 178 +C & .or. itype(i-1).eq.ntyp1 + & ) goto 178 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -3851,6 +3922,7 @@ C & gnmr1(vbld(i),-1.0d0,distchainmax) C else if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then diff = vbld(i)-vbldpDUM +C write(iout,*) i,diff else diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff