'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
#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
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)
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
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
! 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
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)
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. &
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) &
! & +(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
! 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. &
! 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. &
! 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. &
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
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
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
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)
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
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
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) &
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
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
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)
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
!
! 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
! print *,"EVDW KURW",evdw,nres
!!! return
17 continue
+! go to 23
do i=ibond_start,ibond_end
! print *,"I am in EVDW",i
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
+ (( 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,&
!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
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
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)
! 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)
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
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)
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)
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)
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)
projects / unres4.git / blobdiff