if (nfgtasks.gt.1) then
call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
endif
+ if (nres_molec(1).gt.0) then
if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
! write (iout,*) "after make_SCp_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! write (iout,*) "after make_SCSC_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
+ endif
! write (iout,*) "after make_pp_inter_list"
! print *,'Processor',myrank,' calling etotal ipot=',ipot
dGCLdOM1=0.0d0
dPOLdOM1=0.0d0
! write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j
-
+ if (nres_molec(1).eq.0) return
do icont=g_listscsc_start,g_listscsc_end
i=newcontlisti(icont)
j=newcontlistj(icont)
'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
eel_loc=0.0d0
eello_turn3=0.0d0
eello_turn4=0.0d0
+ if (nres_molec(1).eq.0) return
!
if (icheckgrad.eq.1) then
#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
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
! do i=iatscp_s,iatscp_e
+ if (nres_molec(1).eq.0) return
do icont=g_listscp_start,g_listscp_end
i=newcontlistscpi(icont)
j=newcontlistscpj(icont)
iabs(itype(jjj,1)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
-!d write (iout,*) "eij",eij
+! write (iout,*) "eij",eij,iii,jjj
endif
else if (ii.gt.nres .and. jj.gt.nres) then
!c Restraints from contact prediction
itypj=iabs(itype(j,1))
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(nres+j)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) &
+akct*deltad*deltat12 &
+v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr
-! write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
-! & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
-! & " deltat12",deltat12," eij",eij
+! write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, &
+! " akct",akct," deltad",deltad," deltat",deltat1,deltat2, &
+! " deltat12",deltat12," eij",eij
ed=2*akcm*deltad+akct*deltat12
pom1=akct*deltad
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
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
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
#endif
!#define DEBUG
!el write (iout,*) "After sum_gradient"
-#ifdef DEBUG
- write (iout,*) "After sum_gradient"
- do i=1,nres-1
- write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
- write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
- enddo
-#endif
+!#ifdef DEBUG
+! write (iout,*) "After sum_gradient"
+! do i=1,nres-1
+! write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
+! write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
+! enddo
+!#endif
!#undef DEBUG
! If performing constraint dynamics, add the gradients of the constraint energy
if(usampl.and.totT.gt.eq_time) then
! 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
!(3,3,2,maxres)
! allocateion of lists JPRDLA
allocate(newcontlistppi(300*nres))
- allocate(newcontlistscpi(300*nres))
+ allocate(newcontlistscpi(350*nres))
allocate(newcontlisti(300*nres))
allocate(newcontlistppj(300*nres))
- allocate(newcontlistscpj(300*nres))
+ allocate(newcontlistscpj(350*nres))
allocate(newcontlistj(300*nres))
return
real(kind=8) :: facd4, adler, Fgb, facd3
integer troll,jj,istate
real (kind=8) :: dcosom1(3),dcosom2(3)
+ real(kind=8) ::locbox(3)
+ locbox(1)=boxxsize
+ locbox(2)=boxysize
+ locbox(3)=boxzsize
evdw=0.0D0
if (nres_molec(5).eq.0) return
zj=c(3,j)
call to_box(xj,yj,zj)
+! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,xj,yj,zj
+
! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
! aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
! +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
+! write(iout,*) "xj,yj,zj", xj,yj,zj,boxxsize
! dxj = dc_norm( 1, nres+j )
! dyj = dc_norm( 2, nres+j )
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
ctail(k,1)=c(k,i+nres)
ctail(k,2)=c(k,j)
END DO
+ call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+ call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
- Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
- Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
- Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
Rtail = dsqrt( &
(Rtail_distance(1)*Rtail_distance(1)) &
+ (Rtail_distance(2)*Rtail_distance(2)) &
! see unres publications for very informative images
chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
chead(k,2) = c(k, j)
+ 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
! distance
! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
-! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
- Rhead_distance(k) = chead(k,2) - chead(k,1)
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ do k=1,3
+ Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
END DO
! pitagoras (root of sum of squares)
Rhead = dsqrt( &
dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
Fcav = 0.0d0
+ Fisocav=0.0d0
dFdR = 0.0d0
dCAVdOM1 = 0.0d0
dCAVdOM2 = 0.0d0
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
yj=c(2,j)
zj=c(3,j)
call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
dxj = 0.0d0! dc_norm( 1, nres+j )
ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
ctail(k,2)=c(k,j)
END DO
+ call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+ call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
+
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
- Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
- Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
- Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
Rtail = dsqrt( &
(Rtail_distance(1)*Rtail_distance(1)) &
+ (Rtail_distance(2)*Rtail_distance(2)) &
! see unres publications for very informative images
chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
chead(k,2) = c(k, j)
+ ENDDO
! distance
! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
- Rhead_distance(k) = chead(k,2) - chead(k,1)
+ call to_box(chead(1,1),chead(2,1),chead(3,1))
+ call to_box(chead(1,2),chead(2,2),chead(3,2))
+
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ do k=1,3
+ Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
END DO
+
! pitagoras (root of sum of squares)
Rhead = dsqrt( &
(Rhead_distance(1)*Rhead_distance(1)) &
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
yj=c(2,j)
zj=c(3,j)
call to_box(xj,yj,zj)
+! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,xj,yj,zj
! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
! aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
! +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
-
+! write(iout,*) 'after shift', xj,yj,zj
dist_init=xj**2+yj**2+zj**2
itypi=itype(i,2)
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)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
- if ((zi.gt.bordlipbot) &
- .and.(zi.lt.bordliptop)) then
-!C the energy transfer exist
- if (zi.lt.buflipbot) then
-!C what fraction I am in
- fracinbuf=1.0d0- &
- ((zi-bordlipbot)/lipbufthick)
-!C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
+! endif
! print *, sslipi,ssgradlipi
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
zj=c(3,j+nres)
call to_box(xj,yj,zj)
call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
- write(iout,*) "KRUWA", i,j
+! write(iout,*) "KRUWA", i,j
aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
! 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)
include 'mpif.h'
real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
real*8 :: dist_init, dist_temp,r_buff_list
- integer:: contlistscpi(250*nres),contlistscpj(250*nres)
+ integer:: contlistscpi(350*nres),contlistscpj(350*nres)
! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr
projects / unres4.git / blobdiff