! include 'COMMON.TIME1'
real(kind=8) :: time00
!el local variables
- integer :: n_corr,n_corr1,ierror
+ integer :: n_corr,n_corr1,ierror,imatupdate
real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
integer ishield_listbuf(-1:nres), &
shield_listbuf(maxcontsshi,-1:nres),k,j,i,iii,impishi,mojint,jjj
-
-
+! print *,"I START ENERGY"
+ imatupdate=100
+! if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! real(kind=8), dimension(:),allocatable:: fac_shieldbuf
! real(kind=8), dimension(:,:,:),allocatable:: &
! grad_shield_locbuf,grad_shield_sidebuf
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
! call chainbuild_cart
endif
+! print *,"itime_mat",itime_mat,imatupdate
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+ endif
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
+
! print *,'Processor',myrank,' calling etotal ipot=',ipot
! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
etors_nucl=0.0d0
estr_nucl=0.0d0
ecorr3_nucl=0.0d0
+ ecorr_nucl=0.0d0
ebe_nucl=0.0d0
evdwsb=0.0d0
eelsb=0.0d0
eelpsb=0.0d0
evdwpp=0.0d0
eespp=0.0d0
+ etors_d_nucl=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
! print *,"before ecatcat",wcatcat
+ if (nres_molec(5).gt.0) then
if (nfgtasks.gt.1) then
if (fg_rank.eq.0) then
call ecatcat(ecationcation)
else
call ecatcat(ecationcation)
endif
+ if (oldion.gt.0) then
call ecat_prot(ecation_prot)
- call ecats_prot_amber(ecations_prot_amber)
- if (nres_molec(2).gt.0) then
+ else
+ call ecats_prot_amber(ecation_prot)
+ endif
+ else
+ ecationcation=0.0d0
+ ecation_prot=0.0d0
+ endif
+ if ((nres_molec(2).gt.0).and.(nres_molec(1).gt.0)) then
call eprot_sc_base(escbase)
call epep_sc_base(epepbase)
call eprot_sc_phosphate(escpho)
energia(47)=epepbase
energia(48)=escpho
energia(49)=epeppho
- energia(50)=ecations_prot_amber
+! energia(50)=ecations_prot_amber
call sum_energy(energia,.true.)
if (dyn_ss) call dyn_set_nss
! print *," Processor",myrank," left SUM_ENERGY"
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
- ecations_prot_amber=energia(50)
+! ecations_prot_amber=energia(50)
! energia(41)=ecation_prot
! energia(42)=ecationcation
+wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
- +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
+wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
- +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
#endif
energia(0)=etot
! detecting NaNQ
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
- ecations_prot_amber=energia(50)
+! ecations_prot_amber=energia(50)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
estr,wbond,ebe,wang,&
etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
- ecations_prot_amber,etot
+ etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
- ecations_prot_amber,etot
+ etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
! include 'COMMON.SBRIDGE'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj,subchap
+ integer :: iint,itypi,itypi1,itypj,subchap,icont
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dPOLdOM1=0.0d0
- do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
+
+! do i=iatsc_s,iatsc_e
!C print *,"I am in EVDW",i
itypi=iabs(itype(i,1))
! if (i.ne.47) cycle
!
! Calculate SC interaction energy.
!
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+! do iint=1,nint_gr(i)
+! do j=istart(i,iint),iend(i,iint)
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
call dyn_ssbond_ene(i,j,evdwij)
evdw=evdw+evdwij
! Calculate angular part of the gradient.
call sc_grad
ENDIF ! dyn_ss
- enddo ! j
- enddo ! iint
+! enddo ! j
+! enddo ! iint
enddo ! i
! print *,"ZALAMKA", evdw
! write (iout,*) "Number of loop steps in EGB:",ind
#endif
#else
if (i.gt. nnt+2 .and. i.lt.nct+2) then
-! write(iout,*) "i,",molnum(i)
+! write(iout,*) "i,",molnum(i),nloctyp
! print *, "i,",molnum(i),i,itype(i-2,1)
if (molnum(i).eq.1) then
+ if (itype(i-2,1).eq.ntyp1) then
+ iti=nloctyp
+ else
iti = itype2loc(itype(i-2,1))
+ endif
else
iti=nloctyp
endif
0.0d0,1.0d0,0.0d0,&
0.0d0,0.0d0,1.0d0/),shape(unmat))
!el local variables
- integer :: i,k,j
+ integer :: i,k,j,icont
real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
real(kind=8) :: fac,t_eelecij,fracinbuf
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
! print *,"iatel_s,iatel_e,",iatel_s,iatel_e
- do i=iatel_s,iatel_e
+! do i=iatel_s,iatel_e
+! JPRDLC
+ do icont=g_listpp_start,g_listpp_end
+ i=newcontlistppi(icont)
+ j=newcontlistppj(icont)
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
- do j=ielstart(i),ielend(i)
+! do j=ielstart(i),ielend(i)
! write (iout,*) i,j,itype(i,1),itype(j,1)
if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
- enddo ! j
+! enddo ! j
num_cont_hb(i)=num_conti
enddo ! i
! write (iout,*) "Number of loop steps in EELEC:",ind
! sss_ele_grad=0.0d0
! print *,sss_ele_cut,sss_ele_grad,&
! (rij),r_cut_ele,rlamb_ele
-! if (sss_ele_cut.le.0.0) go to 128
+ if (sss_ele_cut.le.0.0) go to 128
rmij=1.0D0/rij
r3ij=rrmij*rmij
!grad enddo
!grad enddo
! 9/28/08 AL Gradient compotents will be summed only at the end
- ggg(1)=facvdw*xj &
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
- ggg(2)=facvdw*yj &
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
- ggg(3)=facvdw*zj &
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj &
*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
do k=1,3
+a32*gmuij1(3)&
+a33*gmuij1(4))&
*fac_shield(i)*fac_shield(j)&
- *sss_ele_cut
+ *sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
!c write(iout,*) "derivative over thatai"
!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
geel_loc_ij*wel_loc&
*fac_shield(i)*fac_shield(j)&
- *sss_ele_cut
+ *sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
!c Derivative over j residue
gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
geel_loc_ji*wel_loc&
*fac_shield(i)*fac_shield(j)&
- *sss_ele_cut
+ *sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
geel_loc_ji=&
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
geel_loc_ji*wel_loc&
*fac_shield(i)*fac_shield(j)&
- *sss_ele_cut
+ *sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
#endif
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
*sss_ele_cut &
*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
*sss_ele_cut &
*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
gacontp_hb1(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
+ ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
- *sss_ele_cut*fac_shield(i)*fac_shield(j)
+ *sss_ele_cut*fac_shield(i)*fac_shield(j) ! &
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gacontp_hb2(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
+ ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
- *sss_ele_cut*fac_shield(i)*fac_shield(j)
+ *sss_ele_cut*fac_shield(i)*fac_shield(j)! &
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gacontp_hb3(k,num_conti,i)=gggp(k) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb1(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
+ ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb2(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
+ ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gacontm_hb3(k,num_conti,i)=gggm(k) &
*sss_ele_cut*fac_shield(i)*fac_shield(j)
+! *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
! Diagnostics. Comment out or remove after debugging!
!C Derivatives in theta
gloc(nphi+i,icg)=gloc(nphi+i,icg) &
+0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3&
- *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)&
+0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
!C#endif
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj,subchap
+ integer :: i,iint,j,k,iteli,itypj,subchap,icont
real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
e1,e2,evdwij,rij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
evdw2_14=0.0d0
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
- do i=iatscp_s,iatscp_e
+! do i=iatscp_s,iatscp_e
+ do icont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(icont)
+ j=newcontlistscpj(icont)
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
- do iint=1,nscp_gr(i)
+! do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+! do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- enddo
+! enddo
- enddo ! iint
+! enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
! & dscp1,dscp2,sumene
! sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
+ if (energy_dec) write (2,*) "i",i," itype",itype(i,1)," it",it, &
+ " escloc",sumene,escloc,it,itype(i,1)
! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
! & ,zz,xx,yy
!#define DEBUG
enddo
#endif
!#undef DEBUG
- do i=1,nres
+ do i=0,nres
do j=1,3
gloc_scbuf(j,i)=gloc_sc(j,i,icg)
enddo
call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
- call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,&
+ call MPI_Reduce(gloc_scbuf(1,0),gloc_sc(1,0,icg),3*nres+3,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
!#define DEBUG
! print *,"gradbuf",gradbufc(1,1),gradc(1,1,icg)
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
- do i=1,nres
+ do i=0,nres
do j=1,1
write (iout,*) i,j,gloc_sc(j,i,icg)
enddo
end subroutine sc_grad
subroutine sc_grad_cat
-! implicit real*8 (a-h,o-z)
use calc_data
-! include 'DIMENSIONS'
-! include 'COMMON.CHAIN'
-! include 'COMMON.DERIV'
-! include 'COMMON.CALC'
-! include 'COMMON.IOUNITS'
real(kind=8), dimension(3) :: dcosom1,dcosom2
-! print *,"wchodze"
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+dCAVdOM1+ dGCLdOM1+ dPOLdOM1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
! eom2=0.0d0
! eom12=evdwij*eps1_om12
! end diagnostics
-! write (iout,*) "eps2der",eps2der," eps3der",eps3der,&
-! " sigder",sigder
-! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
-! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
-!C print *,sss_ele_cut,'in sc_grad'
do k=1,3
dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)&
+ gradpepcatx(k,i)=gradpepcatx(k,i)-gg(k) &
+(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)&
- +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
- +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
+! gradpepcatx(k,j)=gradpepcatx(k,j)+gg(k) &
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
!
! Calculate the components of the gradient in DC and X
!
-!grad do k=i,j-1
-!grad do l=1,3
-!grad gvdwc(l,k)=gvdwc(l,k)+gg(l)
-!grad enddo
-!grad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gradpepcat(l,i)=gradpepcat(l,i)-gg(l)
+ gradpepcat(l,j)=gradpepcat(l,j)+gg(l)
enddo
end subroutine sc_grad_cat
+ subroutine sc_grad_cat_pep
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+ +dCAVdOM1+ dGCLdOM1+ dPOLdOM1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
+ +dCAVdOM2+ dGCLdOM2+ dPOLdOM2
+
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12&
+ +dCAVdOM12+ dGCLdOM12
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+
+ do k=1,3
+ dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) &
+ + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+ *dsci_inv*2.0 &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) &
+ - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) &
+ *dsci_inv*2.0 &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gradpepcat(k,j)=gradpepcat(k,j)+gg(k)
+ enddo
+ end subroutine sc_grad_cat_pep
#ifdef CRYST_THETA
!-----------------------------------------------------------------------------
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-4
+ aincr=1.0D-5
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-7
+ aincr=1.0D-6
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
!#define DEBUG
!el write (iout,*) "After sum_gradient"
#ifdef DEBUG
-!el write (iout,*) "After sum_gradient"
+ 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)
dphi(j,1,i)=0.0d0
dphi(j,2,i)=0.0d0
dphi(j,3,i)=0.0d0
+ dcosomicron(j,1,1,i)=0.0d0
+ dcosomicron(j,1,2,i)=0.0d0
+ dcosomicron(j,2,1,i)=0.0d0
+ dcosomicron(j,2,2,i)=0.0d0
enddo
enddo
! Derivatives of theta's
#else
do i=3,nres
#endif
- if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1).and.molnum(i).ne.5) then
cost1=dcos(omicron(1,i))
sint1=sqrt(1-cost1*cost1)
cost2=dcos(omicron(2,i))
allocate(uygrad(3,3,2,nres))
allocate(uzgrad(3,3,2,nres))
!(3,3,2,maxres)
+! allocateion of lists JPRDLA
+ allocate(newcontlistppi(200*nres))
+ allocate(newcontlistscpi(200*nres))
+ allocate(newcontlisti(200*nres))
+ allocate(newcontlistppj(200*nres))
+ allocate(newcontlistscpj(200*nres))
+ allocate(newcontlistj(200*nres))
return
end subroutine alloc_ener_arrays
!c------------------------------------------------------------------------------
#endif
subroutine ecatcat(ecationcation)
- integer :: i,j,itmp,xshift,yshift,zshift,subchap,k
+ integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj
real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
r7,r4,ecationcation,k0,rcal
real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
epscalc=0.05
r06 = rcat0**6
r012 = r06**2
- k0 = 332.0*(2.0*2.0)/80.0
+! k0 = 332.0*(2.0*2.0)/80.0
itmp=0
do i=1,4
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
-
+! write (iout,*) i,"TUTUT",c(1,i)
+ itypi=itype(i,5)
xi=mod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
yi=mod(yi,boxysize)
if (zi.lt.0) zi=zi+boxzsize
do j=i+1,itmp+nres_molec(5)
+ itypj=itype(j,5)
+! print *,i,j,itypi,itypj
+ k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0
! print *,i,j,'catcat'
xj=c(1,j)
yj=c(2,j)
! r06 = rcat0**6
! r012 = r06**2
! k0 = 332*(2*2)/80
- Evan1cat=epscalc*(r012/rcal**6)
- Evan2cat=epscalc*2*(r06/rcal**3)
+ Evan1cat=epscalc*(r012/(rcal**6))
+ Evan2cat=epscalc*2*(r06/(rcal**3))
Eeleccat=k0/ract
r7 = rcal**7
r4 = rcal**4
gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
enddo
-
+ if (energy_dec) write (iout,*) i,j,Evan1cat,Evan2cat,Eeleccat,&
+ r012,rcal**6,ichargecat(itypi)*ichargecat(itypj)
! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
enddo
end subroutine ecatcat
!---------------------------------------------------------------------------
! new for K+
- subroutine ecats_prot_amber(ecations_prot_amber)
+ subroutine ecats_prot_amber(evdw)
! subroutine ecat_prot2(ecation_prot)
use calc_data
use comm_momo
integer troll,jj,istate
real (kind=8) :: dcosom1(3),dcosom2(3)
- ecations_prot_amber=0.0D0
+ evdw=0.0D0
if (nres_molec(5).eq.0) return
eps_out=80.0d0
! sss_ele_cut=1.0d0
do i=1,4
itmp=itmp+nres_molec(i)
enddo
+! go to 17
! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
do i=ibond_start,ibond_end
! print *,"I am in EVDW",i
itypi=iabs(itype(i,1))
+
! if (i.ne.47) cycle
- if (itypi.eq.ntyp1) cycle
+ if ((itypi.eq.ntyp1).or.(itypi.eq.10)) cycle
itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
do j=itmp+1,itmp+nres_molec(5)
! Calculate SC interaction energy.
- itypj=iabs(itype(j,1))
+ itypj=iabs(itype(j,5))
if ((itypj.eq.ntyp1)) cycle
CALL elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
- dscj_inv=vbld_inv(j+nres)
+ dscj_inv=0.0
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
! alf2 = 0.0d0
! alf12 = 0.0d0
! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
- chi1 = chicat(itypi,itypj)
- chis1 = chiscat(itypi,itypj)
- chip1 = chippcat(itypi,itypj)
+ chi1 = chi1cat(itypi,itypj)
+ chis1 = chis1cat(itypi,itypj)
+ chip1 = chipp1cat(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
! chis2 = chis(itypj,itypi)
-! chis12 = chis1 * chis2
+ chis12 = chis1 * chis2
sig1 = sigmap1cat(itypi,itypj)
! sig2 = sigmap2(itypi,itypj)
! alpha factors from Fcav/Gcav
sigder = -sig * sigsq
rij_shift = 1.0D0 / rij_shift
fac = rij_shift**expon
- c1 = fac * fac * aa_aq(itypi,itypj)
+ c1 = fac * fac * aa_aq_cat(itypi,itypj)
! print *,"ADAM",aa_aq(itypi,itypj)
! c1 = 0.0d0
- c2 = fac * bb_aq(itypi,itypj)
+ c2 = fac * bb_aq_cat(itypi,itypj)
! c2 = 0.0d0
evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
eps2der = eps3rt * evdwij
END DO
erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
erdxj = scalar( ertail(1), dC_norm(1,j) )
- facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
- facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
DO k = 1, 3
pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
- gvdwx(k,i) = gvdwx(k,i) &
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
- (( dFdR + gg(k) ) * pom)
pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
! gvdwx(k,j) = gvdwx(k,j) &
! + (( dFdR + gg(k) ) * pom)
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) &
- (( dFdR + gg(k) ) * ertail(k))
- gvdwc(k,j) = gvdwc(k,j) &
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ (( dFdR + gg(k) ) * ertail(k))
gg(k) = 0.0d0
-
+ ENDDO
!c! Compute head-head and head-tail energies for each state
- isel = iabs(Qi) + iabs(Qj)
+ isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj)
IF (isel.eq.0) THEN
!c! No charges - do nothing
eheadtail = 0.0d0
- ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
+ ELSE IF (isel.eq.1) THEN
!c! Nonpolar-charge interactions
if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
Qi=Qi*2
eheadtail = epol
! eheadtail = 0.0d0
- ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
+ ELSE IF (isel.eq.3) THEN
!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
Qj=Qj*2
Qij=Qij*2
endif
+ write(iout,*) "KURWA0",d1
+
CALL edq_cat(ecl, elj, epol)
eheadtail = ECL + elj + epol
! eheadtail = 0.0d0
- ELSE IF ((isel.eq.2.and. &
- iabs(Qi).eq.1).and. &
- nstate(itypi,itypj).eq.1) THEN
+ ELSE IF ((isel.eq.2)) THEN
!c! Same charge-charge interaction ( +/+ or -/- )
if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
!
! 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
- evdw = evdw + Fcav + eheadtail
+ evdw = evdw + Fcav + eheadtail
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,&
!c!-------------------------------------------------------------------
!c! NAPISY KONCOWE
END DO ! j
- END DO ! iint
END DO ! i
!c write (iout,*) "Number of loop steps in EGB:",ind
!c energy_dec=.false.
! print *,"EVDW KURW",evdw,nres
-
- return
- end subroutine ecats_prot_amber
-
-!---------------------------------------------------------------------------
-! old for Ca2+
- subroutine ecat_prot(ecation_prot)
-! use calc_data
-! use comm_momo
- integer i,j,k,subchap,itmp,inum
- real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
- r7,r4,ecationcation
- real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
- dist_init,dist_temp,ecation_prot,rcal,rocal, &
- Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
- catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
- wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
- costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
- Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
- rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
- opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
- opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
- Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
- ndiv,ndivi
- real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
- gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
- dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
- tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
- v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
- dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
- dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
- dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
- dEvan1Cat
- real(kind=8),dimension(6) :: vcatprm
- ecation_prot=0.0d0
-! first lets calculate interaction with peptide groups
- if (nres_molec(5).eq.0) return
- itmp=0
- do i=1,4
- itmp=itmp+nres_molec(i)
- enddo
-! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+!!! return
+ 17 continue
do i=ibond_start,ibond_end
-! cycle
- if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
- xi=0.5d0*(c(1,i)+c(1,i+1))
- yi=0.5d0*(c(2,i)+c(2,i+1))
- zi=0.5d0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
+
+! print *,"I am in EVDW",i
+ itypi=10 ! the peptide group parameters are for glicine
+
+! if (i.ne.47) cycle
+ if ((itype(i,1).eq.ntyp1).or.itype(i+1,1).eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=dmod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
+ yi=dmod(yi,boxysize)
if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
+ zi=dmod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
-
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=dc_norm(3,i)
+ dsci_inv=vbld_inv(i+1)/2.0
do j=itmp+1,itmp+nres_molec(5)
-! print *,"WTF",itmp,j,i
-! all parameters were for Ca2+ to approximate single charge divide by two
- ndiv=1.0
- if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
- wconst=78*ndiv
- wdip =1.092777950857032D2
- wdip=wdip/wconst
- wmodquad=-2.174122713004870D4
- wmodquad=wmodquad/wconst
- wquad1 = 3.901232068562804D1
- wquad1=wquad1/wconst
- wquad2 = 3
- wquad2=wquad2/wconst
- wvan1 = 0.1
- wvan2 = 6
-! itmp=0
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,5))
+ if ((itypj.eq.ntyp1)) cycle
+ CALL elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=dmod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
xj=xj_temp-xi
yj=yj_temp-yi
zj=zj_temp-zi
- else
+ else
xj=xj_safe-xi
yj=yj_safe-yi
zj=zj_safe-zi
- endif
+ endif
+
+ dxj = 0.0d0! dc_norm( 1, nres+j )
+ dyj = 0.0d0!dc_norm( 2, nres+j )
+ dzj = 0.0d0! dc_norm( 3, nres+j )
+
+ itypi = 10
+ itypj = itype(j,5)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chi1cat(itypi,itypj)
+ chis1 = chis1cat(itypi,itypj)
+ chip1 = chipp1cat(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1cat(itypi,itypj)
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+
+! 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
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! Rtail = 0.0d0
+
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
+ ctail(k,2)=c(k,j)
+ END DO
+!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)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail location and distance calculations
+! dhead1
+ d1 = dheadcat(1, 1, itypi, itypj)
+! print *,"d1",d1
+! d1=0.0d0
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! 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)
+! 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)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_cat(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_cat(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+! print *,"TUT2",fac,chis1,sqom1,pom
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ Chif = Rtail * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+ ENDDO
+!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
+
+ evdw = evdw + Fcav + eheadtail
+
+ 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
+
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_cat_pep
+! END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+
+ return
+ end subroutine ecats_prot_amber
+
+!---------------------------------------------------------------------------
+! old for Ca2+
+ subroutine ecat_prot(ecation_prot)
+! use calc_data
+! use comm_momo
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
! enddo
! enddo
rcpm = sqrt(xj**2+yj**2+zj**2)
real(kind=8) :: facd4, adler, Fgb, facd3
integer troll,jj,istate
real (kind=8) :: dcosom1(3),dcosom2(3)
+ evdw=0.0d0
eps_out=80.0d0
sss_ele_cut=1.0d0
! print *,"EVDW KURW",evdw,nres
facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
- gvdwx(k,i) = gvdwx(k,i) &
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
- dGCLdR * pom&
- dGGBdR * pom&
- dGCVdR * pom&
- dGLJdR * pom
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
- gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
- + dGGBdR * pom+ dGCVdR * pom&
- + dPOLdR1 * (erhead_tail(k,1)&
- -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
- + dPOLdR2 * condor + dGLJdR * pom
+! gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
+! + dGGBdR * pom+ dGCVdR * pom&
+! + dPOLdR1 * (erhead_tail(k,1)&
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
+! + dPOLdR2 * condor + dGLJdR * pom
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) &
- dGCLdR * erhead(k)&
- dGGBdR * erhead(k)&
- dGCVdR * erhead(k)&
- dPOLdR2 * erhead_tail(k,2)&
- dGLJdR * erhead(k)
- gvdwc(k,j) = gvdwc(k,j) &
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ dGCLdR * erhead(k) &
+ dGGBdR * erhead(k) &
+ dGCVdR * erhead(k) &
condor = (erhead_tail(k,2) &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
- gvdwx(k,i) = gvdwx(k,i) &
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
- dPOLdR2 * (erhead_tail(k,2) &
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
- gvdwx(k,j) = gvdwx(k,j) &
- + dPOLdR2 * condor
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dPOLdR2 * condor
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) &
- dPOLdR2 * erhead_tail(k,2)
- gvdwc(k,j) = gvdwc(k,j) &
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ dPOLdR2 * erhead_tail(k,2)
END DO
!c!-------------------------------------------------------------------
!c! ecl
- sparrow = w1 * Qi * om1
- hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
ECL = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ 4.0d0 * hawk / Rhead**5.0d0
!c! dF/dom1
- dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
- dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
!c--------------------------------------------------------------------
!c Polarization energy
!c Epol
!c!-------------------------------------------------------------------
!c! ecl
- sparrow = w1 * Qi * om1
- hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ write(iout,*) "KURWA2",Rhead
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
ECL = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ 4.0d0 * hawk / Rhead**5.0d0
!c! dF/dom1
- dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
- dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
!c--------------------------------------------------------------------
!c--------------------------------------------------------------------
!c Polarization energy
adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
facd1 = d1 * vbld_inv(i+nres)
facd2 = d2 * vbld_inv(j)
- facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
DO k = 1, 3
condor = (erhead_tail(k,2) &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
- gvdwx(k,i) = gvdwx(k,i) &
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
- dGCLdR * pom &
- dPOLdR2 * (erhead_tail(k,2) &
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
- dGLJdR * pom
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
- gvdwx(k,j) = gvdwx(k,j) &
- + dGCLdR * pom &
- + dPOLdR2 * condor &
- + dGLJdR * pom
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
- gvdwc(k,i) = gvdwc(k,i) &
+ gradpepcat(k,i) = gradpepcat(k,i) &
- dGCLdR * erhead(k) &
- dPOLdR2 * erhead_tail(k,2) &
- dGLJdR * erhead(k)
- gvdwc(k,j) = gvdwc(k,j) &
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ dGCLdR * erhead(k) &
+ dPOLdR2 * erhead_tail(k,2) &
+ dGLJdR * erhead(k)
RETURN
END SUBROUTINE edq_cat
+ SUBROUTINE edq_cat_pep(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolcat(itypj,itypi)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+! print *,"CO2", itypi,itypj
+! print *,"CO?!.", w1,w2,Qj,om1
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i) )
+ facd1 = d1 * vbld_inv(i+1)/2.0
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+1)/2.0
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - dGCLdR * pom &
+! - dPOLdR2 * (erhead_tail(k,2) &
+! -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+! - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_cat_pep
SUBROUTINE edd(ECL)
! IMPLICIT NONE
b4cav = alphasur(4,itypi,itypj)
wqd = wquad(itypi, itypj)
!c! used by Fgb
- eps_in = epsintab(itypi,itypj)
+ eps_in = epsintab(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!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)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! 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+nres) + d2 * dc_norm(k, j+nres)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init
+
+
+ SUBROUTINE elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = itype(i,1)
+ itypj = itype(j,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chi1cat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1cat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ 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 )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = ichargecat(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1cat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
-!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
!c!-------------------------------------------------------------------
!c! tail location and distance calculations
Rtail = 0.0d0
DO k = 1, 3
- ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
- ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
END DO
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
!c! Calculate location and distance between polar heads
!c! distance between heads
!c! for each one of our three dimensional space...
- d1 = dhead(1, 1, itypi, itypj)
- d2 = dhead(2, 1, itypi, itypj)
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
DO k = 1,3
!c! location of polar head is computed by taking hydrophobic centre
!c! and moving by a d1 * dc_norm vector
!c! 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+nres) + d2 * dc_norm(k, j+nres)
+ chead(k,2) = c(k, j)
!c! distance
!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
RETURN
- END SUBROUTINE elgrad_init
-
+ END SUBROUTINE elgrad_init_cat
- SUBROUTINE elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ SUBROUTINE elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
use comm_momo
use calc_data
real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
eps_out=80.0d0
- itypi = itype(i,1)
- itypj = itype(j,1)
+ itypi = 10
+ itypj = itype(j,5)
!c! 1/(Gas Constant * Thermostate temperature) = BetaT
!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
!c! t_bath = 300
BetaT = 1.0d0 / (298.0d0 * Rb)
!c! Gay-berne var's
sig0ij = sigmacat( itypi,itypj )
- chi1 = chicat( itypi, itypj )
-! chi2 = chi( itypj, itypi )
+ chi1 = chi1cat( itypi, itypj )
chi2 = 0.0d0
-! chi12 = chi1 * chi2
chi12 = 0.0d0
- chip1 = chippcat( itypi, itypj )
-! chip2 = chipp( itypj, itypi )
+ chip1 = chipp1cat( itypi, itypj )
chip2 = 0.0d0
-! chip12 = chip1 * chip2
chip12 = 0.0d0
-! chi1=0.0
-! chi2=0.0
-! chi12=0.0
-! chip1=0.0
-! chip2=0.0
-! chip12=0.0
!c! not used by momo potential, but needed by sc_angular which is shared
!c! by all energy_potential subroutines
alf1 = 0.0d0
alf2 = 0.0d0
alf12 = 0.0d0
-!c! location, location, location
-! xj = c( 1, nres+j ) - xi
-! yj = c( 2, nres+j ) - yi
-! zj = c( 3, nres+j ) - zi
- 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)
!c! ai*aj from Fgb
- a12sq = rborncat(itypi,itypj) * rborncat(itypj,itypi)
+ a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
!c! a12sq = a12sq * a12sq
!c! charge of amino acid itypi is...
- Qi = ichargecat(itypi)
+ Qi = 0
Qj = ichargecat(itypj)
- Qij = Qi * Qj
+! Qij = Qi * Qj
!c! chis1,2,12
- chis1 = chiscat(itypi,itypj)
-! chis2 = chis(itypj,itypi)
+ chis1 = chis1cat(itypi,itypj)
chis2 = 0.0d0
-! chis12 = chis1 * chis2
chis12 = 0.0d0
sig1 = sigmap1cat(itypi,itypj)
sig2 = sigmap2cat(itypi,itypj)
!c! alpha factors from Fcav/Gcav
b1cav = alphasurcat(1,itypi,itypj)
-! b1cav=0.0
b2cav = alphasurcat(2,itypi,itypj)
b3cav = alphasurcat(3,itypi,itypj)
b4cav = alphasurcat(4,itypi,itypj)
!c! tail location and distance calculations
Rtail = 0.0d0
DO k = 1, 3
- ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
- ctail(k,2)=c(k,j+nres)-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailcat(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
END DO
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
!c! location of polar head is computed by taking hydrophobic centre
!c! and moving by a d1 * dc_norm vector
!c! 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+nres) + d2 * dc_norm(k, j+nres)
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
!c! distance
!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
RETURN
- END SUBROUTINE elgrad_init_cat
-
+ END SUBROUTINE elgrad_init_cat_pep
double precision function tschebyshev(m,n,x,y)
implicit none
return
end function gradtschebyshev
+ subroutine make_SCSC_inter_list
+ 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:: contlisti(200*nres),contlistj(200*nres)
+! integer :: newcontlisti(200*nres),newcontlistj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_sc,g_ilist_sc
+ integer displ(0:nprocs),i_ilist_sc(0:nprocs),ierr
+! print *,"START make_SC"
+ r_buff_list=5.0
+ ilist_sc=0
+ do i=iatsc_s,iatsc_e
+ itypi=iabs(itype(i,1))
+ if (itypi.eq.ntyp1) cycle
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_sc=ilist_sc+1
+! this can be substituted by cantor and anti-cantor
+ contlisti(ilist_sc)=i
+ contlistj(ilist_sc)=j
+
+ endif
+ enddo
+ enddo
+ enddo
+! call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+! MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! call MPI_Gather(newnss,1,MPI_INTEGER,&
+! i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_sc
+ do i=1,ilist_sc
+ write (iout,*) i,contlisti(i),contlistj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_sc,1,MPI_INTEGER,&
+ i_ilist_sc,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_sc(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlisti,ilist_sc,MPI_INTEGER,&
+ newcontlisti,i_ilist_sc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistj,ilist_sc,MPI_INTEGER,&
+ newcontlistj,i_ilist_sc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlisti,g_ilist_sc,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistj,g_ilist_sc,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_sc=ilist_sc
+
+ do i=1,ilist_sc
+ newcontlisti(i)=contlisti(i)
+ newcontlistj(i)=contlistj(i)
+ enddo
+ endif
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_sc
+ do i=1,g_ilist_sc
+ write (iout,*) i,newcontlisti(i),newcontlistj(i)
+ enddo
+#endif
+ call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+ return
+ end subroutine make_SCSC_inter_list
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ subroutine make_SCp_inter_list
+ use MD_data, only: itime_mat
+
+ 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(200*nres),contlistscpj(200*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
+! print *,"START make_SC"
+ r_buff_list=5.0
+ ilist_scp=0
+ do i=iatscp_s,iatscp_e
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do iint=1,nscp_gr(i)
+
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+! Uncomment following three lines for SC-p interactions
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+! Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+#ifdef DEBUG
+ ! r_buff_list is a read value for a buffer
+ if ((sqrt(dist_init).le.(r_cut_ele)).and.(ifirstrun.eq.0)) then
+! Here the list is created
+ ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi_f(ilist_scp_first)=i
+ contlistscpj_f(ilist_scp_first)=j
+ endif
+#endif
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_scp=ilist_scp+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi(ilist_scp)=i
+ contlistscpj(ilist_scp)=j
+ endif
+ enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_scp
+ do i=1,ilist_scp
+ write (iout,*) i,contlistscpi(i),contlistscpj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_scp,g_ilist_scp,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_scp,1,MPI_INTEGER,&
+ i_ilist_scp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_scp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistscpi,ilist_scp,MPI_INTEGER,&
+ newcontlistscpi,i_ilist_scp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistscpj,ilist_scp,MPI_INTEGER,&
+ newcontlistscpj,i_ilist_scp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_scp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistscpi,g_ilist_scp,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistscpj,g_ilist_scp,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_scp=ilist_scp
+
+ do i=1,ilist_scp
+ newcontlistscpi(i)=contlistscpi(i)
+ newcontlistscpj(i)=contlistscpj(i)
+ enddo
+ endif
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_scp
+ do i=1,g_ilist_scp
+ write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
+ enddo
+
+! if (ifirstrun.eq.0) ifirstrun=1
+! do i=1,ilist_scp_first
+! do j=1,g_ilist_scp
+! if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+! (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+! enddo
+! print *,itime_mat,"ERROR matrix needs updating"
+! print *,contlistscpi_f(i),contlistscpj_f(i)
+! 126 continue
+! enddo
+#endif
+ call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
+ return
+ end subroutine make_SCp_inter_list
+
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
+
+ subroutine make_pp_inter_list
+ 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 :: xmedj,ymedj,zmedj
+ real*8 :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
+ real*8 :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
+ integer:: contlistppi(200*nres),contlistppj(200*nres)
+! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
+ integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
+! print *,"START make_SC"
+ ilist_pp=0
+ r_buff_list=5.0
+ do i=iatel_s,iatel_e
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ do j=ielstart(i),ielend(i)
+! write (iout,*) i,j,itype(i,1),itype(j,1)
+ if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
+
+! 1,j)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ endif
+ enddo
+ enddo
+ enddo
+
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_pp=ilist_pp+1
+! this can be substituted by cantor and anti-cantor
+ contlistppi(ilist_pp)=i
+ contlistppj(ilist_pp)=j
+ endif
+ enddo
+ enddo
+! enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_pp
+ do i=1,ilist_pp
+ write (iout,*) i,contlistppi(i),contlistppj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_pp,g_ilist_pp,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_pp,1,MPI_INTEGER,&
+ i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_pp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER,&
+ newcontlistppi,i_ilist_pp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER,&
+ newcontlistppj,i_ilist_pp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_pp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM,IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+ g_ilist_pp=ilist_pp
+
+ do i=1,ilist_pp
+ newcontlistppi(i)=contlistppi(i)
+ newcontlistppj(i)=contlistppj(i)
+ enddo
+ endif
+ call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_pp
+ do i=1,g_ilist_pp
+ write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
+ enddo
+#endif
+ return
+ end subroutine make_pp_inter_list
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff