edihcnstr=0
endif
- if (constr_homology.ge.1) then
+ if (constr_homology.ge.1.and.waga_homology(iset).ne.0d0) then
call e_modeller(ehomology_constr)
c print *,'iset=',iset,'me=',me,ehomology_constr,
c & 'Processor',fg_rank,' CG group',kolor,
#endif
enddo
enddo
- if (constr_homology.gt.0) then
+ if (constr_homology.gt.0.and.waga_homology(iset).ne.0d0) then
do i=1,nct
do j=1,3
gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
-
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+ call flush(iout)
+ endif
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
C Cartesian derivatives
+!DIR$ UNROLL(0)
do l=1,3
c ghalf1=0.5d0*agg(l,1)
c ghalf2=0.5d0*agg(l,2)
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
else if (ii.gt.nres .and. jj.gt.nres) then
c Restraints from contact prediction
dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)')
+ & "edisl",ii,jj,
+ & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+ & fordepth(i),dd
+ else
if (dhpb1(i).gt.0.0d0) then
ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
C Evaluate gradient.
C
fac=waga*rdis/dd
- endif
+ endif
+ endif
do j=1,3
ggg(j)=fac*(c(j,jj)-c(j,ii))
enddo
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)')
+ 7 "edisl",ii,jj,
+ & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+ & fordepth(i),dd
+c if (energy_dec)
+c & write (iout,*) fac
+ else
if (dhpb1(i).gt.0.0d0) then
ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
C
fac=waga*rdis/dd
endif
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
do j=1,3
enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+c do i=1,nres
+c write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3)
+c enddo
return
end
C--------------------------------------------------------------------------
do i=ibondp_start,ibondp_end
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
diff=vbld(i+nres)-vbldsc0(1,iti)
c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
c & AKSC(1,iti),AKSC(1,iti)*diff*diff
+ if (energy_dec) then
+ write (iout,*)
+ & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+ & AKSC(1,iti),AKSC(1,iti)*diff*diff
+ call flush(iout)
+ endif
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
sinkt(k)=dsin(k*theti2)
enddo
C if (i.gt.3) then
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
& phii1*rad2deg,ethetai
c lprn1=.false.
etheta=etheta+ethetai
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'ebend',i,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)=gloc(nphi+i-2,icg)+wang*dethetai
c
- do i=1,19
+ do i=1,max_template
distancek(i)=9999999.9
enddo
j = jres_homo(ii)
dij=dist(i,j)
c write (iout,*) "dij(",i,j,") =",dij
+ nexl=0
do k=1,constr_homology
+c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
+ if(.not.l_homo(k,ii)) then
+ nexl=nexl+1
+ cycle
+ endif
distance(k)=odl(k,ii)-dij
c write (iout,*) "distance(",k,") =",distance(k)
c
endif
enddo
- min_odl=minval(distancek)
+
+c min_odl=minval(distancek)
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
+ & min_odl=distancek(kk)
+ enddo
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
write (iout,* )"min_odl",min_odl
#endif
+#ifdef OLDRESTR
odleg2=0.0d0
+#else
+ if (waga_dist.ge.0.0d0) then
+ odleg2=nexl
+ else
+ odleg2=0.0d0
+ endif
+#endif
do k=1,constr_homology
c Nie wiem po co to liczycie jeszcze raz!
c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
c & (2*(sigma_odl(i,j,k))**2))
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c
c For Gaussian-type Urestr
c & *waga_dist)+min_odl
c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
c
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c For Gaussian-type Urestr
c
enddo
#endif
do i=idihconstr_start_homo,idihconstr_end_homo
+#ifdef OLDRESTR
kat2=0.0d0
+#else
+ kat2=nexl
+#endif
c betai=beta(i,i+1,i+2,i+3)
- betai = phi(i+3)
+ betai = phi(i)
c write (iout,*) "betai =",betai
do k=1,constr_homology
dih_diff(k)=pinorm(dih(k,i)-betai)
c & -(6.28318-dih_diff(i,k))
c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
c & 6.28318+dih_diff(i,k)
-
+#ifdef OLD_DIHED
kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+ kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#endif
c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
gdih(k)=dexp(kat3)
kat2=kat2+gdih(k)
sum_gdih=kat2
sum_sgdih=0.0d0
do k=1,constr_homology
+#ifdef OLD_DIHED
sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+#else
+ sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd
+#endif
c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
sum_sgdih=sum_sgdih+sgdih
enddo
c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
ccc & gloc(nphi+i-3,icg)
- gloc(i,icg)=gloc(i,icg)+grad_dih3
+ gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
c if (i.eq.25) then
c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
c endif
c Deviation of theta angles wrt constr_homology ref structures
c
utheta_i=0.0d0 ! argument of Gaussian for single k
+#ifdef OLDRESTR
gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+#else
+ gutheta_i=nexl
+#endif
c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
c over residues in a fragment
c write (iout,*) "theta(",i,")=",theta(i)
utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
- gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk)
+ gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk)
c Gradient for single Gaussian restraint in subr Econstr_back
c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
c
#endif
do i=loc_start,loc_end
usc_diff_i=0.0d0 ! argument of Gaussian for single k
+#ifdef OLDRESTR
guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+#else
+ guscdiff(i)=nexl
+#endif
c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
c write(iout,*) "xxtab, yytab, zztab"
c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
c uscdiffk(k)=usc_diff(i)
guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
- guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk)
+ guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk)
c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
c & xxref(j),yyref(j),zzref(j)
enddo
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
etors_d=0.0D0
do i=iphid_start,iphid_end
+ etors_d_ii=0.0D0
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
sinphi2=dsin(j*phii1)
etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
& v2cij*cosphi2+v2sij*sinphi2
+ if (energy_dec) etors_d_ii=etors_d_ii+
+ & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
sinphi1m2=dsin(l*phii-(k-l)*phii1)
etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
& v1sdij*sinphi1p2+v2sdij*sinphi1m2
+ if (energy_dec) etors_d_ii=etors_d_ii+
+ & v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2
gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
& -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
& -v1cdij*sinphi1p2+v2cdij*sinphi1m2)
enddo
enddo
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'etor_d',i,etors_d_ii
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
c write (iout,*) "gloci", gloc(i-3,icg)
projects / unres.git / blobdiff