X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=ec3ffd409508b68b9cc76e7e12856de4b24c1333;hb=b2ef3fdfb9fd19710ea49377832d6c3791f04c11;hp=58f17148294183bad567aa60a66274276aaaafba;hpb=1d7c7ce298528a7de69c6dced5a52e9780c9ce12;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 58f1714..ec3ffd4 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -1953,7 +1953,7 @@ 'evdw',i,j,evdwij,' ss' ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij,' ss' - do k=j+1,iend(i,iint) + do k=j+1,nres !C search over all next residues if (dyn_ss_mask(k)) then !C check if they are cysteins @@ -5064,14 +5064,19 @@ #ifdef NEWCORR gloc(nphi+i,icg)=gloc(nphi+i,icg)& -(gs13+gsE13+gsEE1)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)& -(gs23+gs21+gsEE2)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)& -(gs32+gsE31+gsEE3)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + !c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- !c & gs2 @@ -11765,7 +11770,7 @@ enddo do k=1,3 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)) -!C print *,'gg',k,gg(k) +! print *,'gg',k,gg(k) enddo ! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut ! write (iout,*) "gg",(gg(k),k=1,3) @@ -12521,6 +12526,12 @@ do j=1,3 grad_s(j,i)=gcart(j,i) grad_s(j+3,i)=gxcart(j,i) + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + enddo enddo else @@ -12745,12 +12756,15 @@ do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) -! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) - -! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + else !- split gradient check call zerograd @@ -16761,8 +16775,8 @@ chip1=chip(itypi) ! if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg) enddo #ifdef DEBUG - write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),& - (gxcart(j,i),j=1,3),gloc(i,icg) + write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3),gloc(i,icg),(gloc_sc(j,i,icg),j=1,3) #endif enddo #ifdef TIMING @@ -17091,10 +17105,12 @@ chip1=chip(itypi) do j=1,3 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/& vbld(i-1) - if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint + if (((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))) & + dtheta(j,1,i)=-dcostheta(j,1,i)/sint dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/& vbld(i) - if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint + if ((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))& + dtheta(j,2,i)=-dcostheta(j,2,i)/sint enddo enddo #if defined(MPI) && defined(PARINTDER) @@ -17151,11 +17167,23 @@ chip1=chip(itypi) cost1=dcos(theta(i-1)) cosg=dcos(phi(i)) scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. & phi(i).gt.pi34.and.phi(i).le.pi.or. & @@ -17164,8 +17192,16 @@ chip1=chip(itypi) call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2) call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) do j=1,3 + if (sint.ne.0.0d0) then ctgt=cost/sint + else + ctgt=0.0d0 + endif + if (sint1.ne.0.0d0) then ctgt1=cost1/sint1 + else + ctgt1=0.0d0 + endif cosg_inv=1.0d0/cosg if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) & @@ -17180,6 +17216,10 @@ chip1=chip(itypi) ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) dphi(j,3,i)=cosg_inv*dsinphi(j,3,i) endif +! write(iout,*) "just after,close to pi",dphi(j,3,i),& +! sing*(ctgt1*dtheta(j,2,i-1)),ctgt*dtheta(j,1,i), & +! (fac0*vp2(j)+sing*dc_norm(j,i-2)),vbld_inv(i-1) + ! Bug fixed 3/24/05 (AL) enddo ! Obtaining the gamma derivatives from cosine derivative @@ -17234,12 +17274,25 @@ chip1=chip(itypi) ! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1)) + ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) - ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4 + else + fac4=0.0d0 + endif + ! Obtaining the gamma derivatives from sine derivative if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. & @@ -17307,11 +17360,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) ! enddo scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or. & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or. & @@ -17382,11 +17447,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. & @@ -19044,14 +19121,12 @@ chip1=chip(itypi) if (idssb(i).eq.newihpb(j) .and. & jdssb(i).eq.newjhpb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i) #endif -#endif enddo do i=1,newnss @@ -19061,21 +19136,22 @@ chip1=chip(itypi) if (newihpb(i).eq.idssb(j) .and. & newjhpb(i).eq.jdssb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) #endif -#endif enddo - +!#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) nss=newnss do i=1,nss idssb(i)=newihpb(i) jdssb(i)=newjhpb(i) enddo +!#else +! nss=0 +!#endif return end subroutine dyn_set_nss @@ -22371,6 +22447,11 @@ chip1=chip(itypi) b3cav = alphasurcat(3,itypi,itypj) b4cav = alphasurcat(4,itypi,itypj) +! b1cav=0.0d0 +! b2cav=0.0d0 +! b3cav=0.0d0 +! b4cav=0.0d0 + ! used to determine whether we want to do quadrupole calculations eps_in = epsintabcat(itypi,itypj) if (eps_in.eq.0.0) eps_in=1.0 @@ -22406,7 +22487,7 @@ chip1=chip(itypi) enddo call to_box(chead(1,1),chead(2,1),chead(3,1)) call to_box(chead(1,2),chead(2,2),chead(3,2)) - write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1 +! write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1 ! distance ! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) ! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) @@ -22459,6 +22540,16 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 + if (evdw.gt.1.0d6) then + write (*,'(2(1x,a3,i3),7f7.2)') & + restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& + 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq + write(*,*) facsig,faceps1_inv,om1,chiom1,chi1 + write(*,*) "ANISO?!",chi1 +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw + endif + RETURN END IF sigder = -sig * sigsq @@ -22492,7 +22583,7 @@ chip1=chip(itypi) gg(1) = fac gg(2) = fac gg(3) = fac - +! print *,"GG(1),distance grad",gg(1) fac = chis1 * sqom1 + chis2 * sqom2 & - 2.0d0 * chis12 * om1 * om2 * om12 pom = 1.0d0 - chis1 * chis2 * sqom12 @@ -22531,12 +22622,12 @@ chip1=chip(itypi) erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) erdxj = scalar( ertail(1), dC_norm(1,j) ) facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) - facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres) + facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j) DO k = 1, 3 pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) gradpepcatx(k,i) = gradpepcatx(k,i) & - (( dFdR + gg(k) ) * pom) - pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) + pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j)) ! gvdwx(k,j) = gvdwx(k,j) & ! + (( dFdR + gg(k) ) * pom) gradpepcat(k,i) = gradpepcat(k,i) & @@ -22546,7 +22637,10 @@ chip1=chip(itypi) gg(k) = 0.0d0 ENDDO !c! Compute head-head and head-tail energies for each state +!! if (.false.) then ! turn off electrostatic + if (itype(j,5).gt.0) then ! the normal cation case isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj) +! print *,i,itype(i,1),isel IF (isel.eq.0) THEN !c! No charges - do nothing eheadtail = 0.0d0 @@ -22557,10 +22651,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL enq_cat(epol) eheadtail = epol @@ -22572,10 +22662,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif ! write(iout,*) "KURWA0",d1 CALL edq_cat(ecl, elj, epol) @@ -22589,10 +22675,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) eheadtail = ECL + Egb + Epol + Fisocav + Elj @@ -22613,16 +22695,28 @@ chip1=chip(itypi) ! ! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav + else + write(iout,*) "not yet implemented",j,itype(j,5) + endif +!! endif ! turn off electrostatic evdw = evdw + Fcav + eheadtail +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& Equad,evdwij+Fcav+eheadtail,evdw ! evdw = evdw + Fcav + eheadtail - +! print *,"before sc_grad_cat", i,j, gradpepcat(1,j) ! iF (nstate(itypi,itypj).eq.1) THEN CALL sc_grad_cat +! print *,"after sc_grad_cat", i,j, gradpepcat(1,j) + ! END IF !c!------------------------------------------------------------------- !c! NAPISY KONCOWE @@ -22633,6 +22727,7 @@ chip1=chip(itypi) ! print *,"EVDW KURW",evdw,nres !!! return 17 continue +! go to 23 do i=ibond_start,ibond_end ! print *,"I am in EVDW",i @@ -22794,6 +22889,13 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),6f6.2)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif RETURN END IF sigder = -sig * sigsq @@ -22885,24 +22987,26 @@ chip1=chip(itypi) + (( dFdR + gg(k) ) * ertail(k)) gg(k) = 0.0d0 ENDDO + if (itype(j,5).gt.0) then !c! Compute head-head and head-tail energies for each state isel = 3 !c! Dipole-charge interactions - if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then - Qi=Qi*2 - Qij=Qij*2 - endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL edq_cat_pep(ecl, elj, epol) eheadtail = ECL + elj + epol ! print *,"i,",i,eheadtail ! eheadtail = 0.0d0 - + else +!HERE WATER and other types of molecules solvents will be added + write(iout,*) "not yet implemented" +! CALL edd_cat_pep + endif evdw = evdw + Fcav + eheadtail - +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& @@ -22919,7 +23023,8 @@ chip1=chip(itypi) !c write (iout,*) "Number of loop steps in EGB:",ind !c energy_dec=.false. ! print *,"EVDW KURW",evdw,nres - + 23 continue +! print *,"before leave sc_grad_cat", i,j, gradpepcat(1,nres-1) return end subroutine ecats_prot_amber @@ -23474,9 +23579,13 @@ chip1=chip(itypi) dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, & dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, & dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, & - dEcavdCm + dEcavdCm,boxik real(kind=8),dimension(14) :: vcatnuclprm ecation_nucl=0.0d0 + boxik(1)=boxxsize + boxik(2)=boxysize + boxik(3)=boxzsize + if (nres_molec(5).eq.0) return itmp=0 do i=1,4 @@ -23519,12 +23628,16 @@ chip1=chip(itypi) vsug(k)=c(k,i) vcat(k)=c(k,j) enddo + call to_box(vcm(1),vcm(2),vcm(3)) + call to_box(vsug(1),vsug(2),vsug(3)) + call to_box(vcat(1),vcat(2),vcat(3)) do k=1,3 - dx(k) = vcat(k)-vcm(k) - enddo - do k=1,3 +! dx(k) = vcat(k)-vcm(k) +! enddo + dx(k)=boxshift(vcat(k)-vcm(k),boxik(k)) +! do k=1,3 v1(k)=dc(k,i+nres) - v2(k)=(vcat(k)-vsug(k)) + v2(k)=boxshift(vcat(k)-vsug(k),boxik(k)) enddo v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2) v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3) @@ -25920,7 +26033,7 @@ chip1=chip(itypi) ! integer :: k !c! Epol and Gpol analytical parameters alphapol1 = alphapolcat(itypi,itypj) - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat2(itypj,itypi) !c! Fisocav and Gisocav analytical parameters al1 = alphisocat(1,itypi,itypj) al2 = alphisocat(2,itypi,itypj) @@ -26573,7 +26686,7 @@ chip1=chip(itypi) use calc_data use comm_momo double precision facd3, adler,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) !c! R2 - distance between head of jth side chain and tail of ith sidechain R2 = 0.0d0 DO k = 1, 3 @@ -26871,7 +26984,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -26990,7 +27103,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -27336,9 +27449,9 @@ chip1=chip(itypi) alf1 = 0.0d0 alf2 = 0.0d0 alf12 = 0.0d0 - dxj = dc_norm( 1, nres+j ) - dyj = dc_norm( 2, nres+j ) - dzj = dc_norm( 3, nres+j ) + dxj = 0.0d0 !dc_norm( 1, nres+j ) + dyj = 0.0d0 !dc_norm( 2, nres+j ) + dzj = 0.0d0 !dc_norm( 3, nres+j ) !c! distance from center of chain(?) to polar/charged head d1 = dheadcat(1, 1, itypi, itypj) d2 = dheadcat(2, 1, itypi, itypj) @@ -27588,8 +27701,10 @@ chip1=chip(itypi) zi=c(3,nres+i) call to_box(xi,yi,zi) do iint=1,nint_gr(i) +! print *,"is it wrong", iint,i do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j,1)) + if (energy_dec) write(iout,*) "LISTA ZAKRES",istart(i,iint),iend(i,iint),iatsc_s,iatsc_e if (itypj.eq.ntyp1) cycle xj=c(1,nres+j) yj=c(2,nres+j)