X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=030de74baa63e988b875036a7d3af0aab6e0c3b7;hb=a7766281a78077166f3389f4dd9e03c98953f6c3;hp=c02d085d8739b55738e1b2f27086aa05f6ebc34d;hpb=51c67c951b6b9f135cda703cc18d60030aff6129;p=unres.git diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index c02d085..030de74 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -363,9 +363,9 @@ cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -378,7 +378,7 @@ C cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -528,9 +528,9 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -539,7 +539,7 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -632,9 +632,9 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -648,7 +648,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -761,9 +761,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -777,7 +777,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -900,9 +900,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=iabs(itype(i)) if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -916,7 +916,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -2769,7 +2769,7 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) if (itypj.eq.21) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi @@ -2880,7 +2880,8 @@ C iii and jjj point to the residues for which the distance is assigned. endif 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 + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij else @@ -2940,7 +2941,7 @@ C include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=itype(i) + itypi=iabs(itype(i)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -2948,7 +2949,7 @@ C dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=dsc_inv(itypi) - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=dsc_inv(itypj) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -3026,6 +3027,7 @@ c logical energy_dec /.false./ double precision u(3),ud(3) estr=0.0d0 + estr1=0.0d0 do i=nnt+1,nct if (itype(i-1).eq.21 .or. itype(i).eq.21) then estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) @@ -3045,12 +3047,12 @@ c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff endif enddo - estr=0.5d0*AKP*estr + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=nnt,nct - iti=itype(i) + iti=iabs(itype(i)) if (iti.ne.10 .and. iti.ne.21) then nbi=nbondterm(iti) if (nbi.eq.1) then @@ -3131,6 +3133,18 @@ c write (iout,*) ithet_start,ithet_end C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) + ichir1=isign(1,itype(i-2)) + ichir2=isign(1,itype(i)) + if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) + if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) + if (itype(i-1).eq.10) then + itype1=isign(10,itype(i-2)) + ichir11=isign(1,itype(i-2)) + ichir12=isign(1,itype(i-2)) + itype2=isign(10,itype(i)) + ichir21=isign(1,itype(i)) + ichir22=isign(1,itype(i)) + endif if (i.gt.3 .and. itype(i-2).ne.21) then #ifdef OSF phii=phi(i) @@ -3168,8 +3182,12 @@ C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2). C In following comments this theta will be referred to as t_c. thet_pred_mean=0.0d0 do k=1,2 - athetk=athet(k,it) - bthetk=bthet(k,it) + athetk=athet(k,it,ichir1,ichir2) + bthetk=bthet(k,it,ichir1,ichir2) + if (it.eq.10) then + athetk=athet(k,itype1,ichir11,ichir12) + bthetk=bthet(k,itype2,ichir21,ichir22) + endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) enddo c write (iout,*) "thet_pred_mean",thet_pred_mean @@ -3177,8 +3195,16 @@ c write (iout,*) "thet_pred_mean",thet_pred_mean thet_pred_mean=thet_pred_mean*ss+a0thet(it) c write (iout,*) "thet_pred_mean",thet_pred_mean C Derivatives of the "mean" values in gamma1 and gamma2. - dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss - dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss + dthetg1=(-athet(1,it,ichir1,ichir2)*y(2) + &+athet(2,it,ichir1,ichir2)*y(1))*ss + dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2) + & +bthet(2,it,ichir1,ichir2)*z(1))*ss + if (it.eq.10) then + dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2) + &+athet(2,itype1,ichir11,ichir12)*y(1))*ss + dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2) + & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss + endif if (theta(i).gt.pi-delta) then call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, & E_tc0) @@ -3354,7 +3380,8 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(itype(i-1)) +CC Ta zmina jest niewlasciwa + ityp2=ithetyp(iabs(itype(i-1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) @@ -3366,7 +3393,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) #else phii=phi(i) #endif - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp(iabs(itype(i-2))) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@ -3387,7 +3414,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) #else phii1=phi(i+1) #endif - ityp3=ithetyp(itype(i)) + ityp3=ithetyp(iabs(itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@ -3542,7 +3569,7 @@ c write (iout,'(a)') 'ESC' it=itype(i) if (it.eq.21) cycle if (it.eq.10) goto 1 - nlobit=nlob(it) + nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol @@ -3697,7 +3724,7 @@ C Compute the contribution to SC energy and derivatives do iii=-1,1 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac do k=1,3 @@ -3778,7 +3805,7 @@ C Compute the contribution to SC energy and derivatives dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@ -3842,7 +3869,7 @@ C cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=iabs(itype(i)) if (it.eq.10) goto 1 c C Compute the axes of tghe local cartesian coordinates system; store in @@ -3882,7 +3909,7 @@ c do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + z_prime(j)*dc_norm(j,i+nres) + zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@ -3892,7 +3919,7 @@ C C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@ -4312,14 +4339,19 @@ c lprn=.true. if (itype(i-2).eq.21 .or. itype(i-1).eq.21 & .or. itype(i).eq.21) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 + if (iabs(itype(i)).eq.20) then + iblock=2 + else + iblock=1 + endif itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1) - v1ij=v1(j,itori,itori1) - v2ij=v2(j,itori,itori1) + do j=1,nterm(itori,itori1,iblock) + v1ij=v1(j,itori,itori1,iblock) + v2ij=v2(j,itori,itori1,iblock) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@ -4332,7 +4364,7 @@ C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1 C cosphi=dcos(0.5d0*phii) sinphi=dsin(0.5d0*phii) - do j=1,nlor(itori,itori1) + do j=1,nlor(itori,itori1,iblock) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@ -4343,11 +4375,11 @@ C gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1) + etors=etors-v0(itori,itori1,iblock) 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, - & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) 1215 continue @@ -4414,12 +4446,14 @@ c lprn=.true. phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 + iblock=1 + if (iabs(itype(i+1)).eq.20) iblock=2 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) - v2cij=v1c(2,j,itori,itori1,itori2) - v2sij=v1s(2,j,itori,itori1,itori2) + do j=1,ntermd_1(itori,itori1,itori2,iblock) + v1cij=v1c(1,j,itori,itori1,itori2,iblock) + v1sij=v1s(1,j,itori,itori1,itori2,iblock) + v2cij=v1c(2,j,itori,itori1,itori2,iblock) + v2sij=v1s(2,j,itori,itori1,itori2,iblock) cosphi1=dcos(j*phii) sinphi1=dsin(j*phii) cosphi2=dcos(j*phii1) @@ -4429,12 +4463,12 @@ C Regular cosine and sine terms gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo - do k=2,ntermd_2(itori,itori1,itori2) + do k=2,ntermd_2(itori,itori1,itori2,iblock) do l=1,k-1 - v1cdij = v2c(k,l,itori,itori1,itori2) - v2cdij = v2c(l,k,itori,itori1,itori2) - v1sdij = v2s(k,l,itori,itori1,itori2) - v2sdij = v2s(l,k,itori,itori1,itori2) + v1cdij = v2c(k,l,itori,itori1,itori2,iblock) + v2cdij = v2c(l,k,itori,itori1,itori2,iblock) + v1sdij = v2s(k,l,itori,itori1,itori2,iblock) + v2sdij = v2s(l,k,itori,itori1,itori2,iblock) cosphi1p2=dcos(l*phii+(k-l)*phii1) cosphi1m2=dcos(l*phii-(k-l)*phii1) sinphi1p2=dsin(l*phii+(k-l)*phii1)