gvdwc_peppho
!------------------------------IONS GRADIENT
real(kind=8),dimension(:,:),allocatable :: gradcatcat, &
- gradpepcat,gradpepcatx
+ gradpepcat,gradpepcatx,gradnuclcat,gradnuclcatx
! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
Eafmforce,ethetacnstr
- real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
+ real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6,ehomology_constr
! now energies for nulceic alone parameters
real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
! energies for ions
- real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+ ecation_nucl
! energies for protein nucleic acid interaction
real(kind=8) :: escbase,epepbase,escpho,epeppho
! allocate(ishield_listbuf(nres))
! allocate(shield_listbuf(maxcontsshi,nres))
! endif
-
+! print *,"wstrain check", wstrain
! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
! & " nfgtasks",nfgtasks
if (nfgtasks.gt.1) then
weights_(47)=wpepbase
weights_(48)=wscpho
weights_(49)=wpeppho
+ weights_(50)=wcatnucl
! wcatcat= weights(41)
! wcatprot=weights(42)
wpepbase=weights(47)
wscpho=weights(48)
wpeppho=weights(49)
+ wcatnucl=weights(50)
! welpsb=weights(28)*fact(1)
!
! wcorr_nucl= weights(37)*fact(1)
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"
+! 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"
+! write (iout,*) "after make_SCSC_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
- write (iout,*) "after make_pp_inter_list"
+ endif
+! write (iout,*) "after make_pp_inter_list"
! print *,'Processor',myrank,' calling etotal ipot=',ipot
! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
!c print *,"Processor",myrank," computed Utor"
! print *,"Processor",myrank," computed Utor"
-
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+! print *,'iset=',iset,'me=',me,ehomology_constr,
+! & 'Processor',fg_rank,' CG group',kolor,
+! & ' absolute rank',MyRank
+! print *,"tu"
+ else
+ ehomology_constr=0.0d0
+ endif
+
!
! 6/23/01 Calculate double-torsional energy
!
!
! If performing constraint dynamics, call the constraint energy
! after the equilibration time
- if(usampl.and.totT.gt.eq_time) then
-!elwrite(iout,*) "afeter multibody hb"
+ if((usampl).and.(totT.gt.eq_time)) then
+ write(iout,*) "usampl",usampl
call EconstrQ
!elwrite(iout,*) "afeter multibody hb"
call Econstr_back
call epsb(evdwpsb,eelpsb)
call esb(esbloc)
call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+ call ecat_nucl(ecation_nucl)
else
etors_nucl=0.0d0
estr_nucl=0.0d0
evdwpp=0.0d0
eespp=0.0d0
etors_d_nucl=0.0d0
+ ecation_nucl=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
! print *,"before ecatcat",wcatcat
energia(48)=escpho
energia(49)=epeppho
! energia(50)=ecations_prot_amber
+ energia(50)=ecation_nucl
+ energia(51)=ehomology_constr
call sum_energy(energia,.true.)
if (dyn_ss) call dyn_set_nss
! print *," Processor",myrank," left SUM_ENERGY"
eliptran,etube, Eafmforce,ethetacnstr
real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
- ecorr3_nucl
- real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+ ecorr3_nucl,ehomology_constr
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+ ecation_nucl
real(kind=8) :: escbase,epepbase,escpho,epeppho
integer :: i
#ifdef MPI
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
+ ecation_nucl=energia(50)
+ ehomology_constr=energia(51)
! ecations_prot_amber=energia(50)
! energia(41)=ecation_prot
+wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
+wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
- +Eafmforce+ethetacnstr &
+ +Eafmforce+ethetacnstr+ehomology_constr &
+wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+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
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
+wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
+wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
- +Eafmforce+ethetacnstr &
+ +Eafmforce+ethetacnstr+ehomology_constr &
+wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+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
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl
#endif
energia(0)=etot
! detecting NaNQ
etube,ethetacnstr,Eafmforce
real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
- ecorr3_nucl
- real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+ ecorr3_nucl,ehomology_constr
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+ ecation_nucl
real(kind=8) :: escbase,epepbase,escpho,epeppho
etot=energia(0)
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
+ ecation_nucl=energia(50)
+ ehomology_constr=energia(51)
+
! ecations_prot_amber=energia(50)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
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,&
- etot
+ ecation_nucl,wcatnucl,ehomology_constr,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
+ 'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/&
+ 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/&
'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce, &
- etube,wtube, &
+ etube,wtube, ehomology_constr,&
estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
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,&
- etot
+ ecation_nucl,wcatnucl,ehomology_constr,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
+ 'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/&
+ 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/&
'ETOT= ',1pE16.6,' (total)')
#endif
return
dCAVdOM1=0.0d0
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)
-
+! write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j
! do i=iatsc_s,iatsc_e
!C print *,"I am in EVDW",i
itypi=iabs(itype(i,1))
'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
zj=c(3,nres+j)
call to_box(xj,yj,zj)
call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! write (iout,*) "KWA2", itypi,itypj
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 &
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
etors_d=0.0d0
return
end subroutine etor_d
+!-----------------------------------------------------------------------------
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+ subroutine e_modeller(ehomology_constr)
+ real(kind=8) :: ehomology_constr
+ ehomology_constr=0.0d0
+ write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!"
+ return
+ end subroutine e_modeller
+C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
#else
!-----------------------------------------------------------------------------
subroutine etor(etors)
return
end subroutine etor_d
#endif
+!----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+ subroutine e_modeller(ehomology_constr)
+! implicit none
+! include 'DIMENSIONS'
+ use MD_data, only: iset
+ real(kind=8) :: ehomology_constr
+ integer nnn,i,ii,j,k,ijk,jik,ki,kk,nexl,irec,l
+ integer katy, odleglosci, test7
+ real(kind=8) :: odleg, odleg2, odleg3, kat, kat2, kat3
+ real(kind=8) :: Eval,Erot,min_odl
+ real(kind=8),dimension(constr_homology) :: distance,distancek,godl,dih_diff,gdih, &
+ gtheta,dscdiff, &
+ uscdiffk,guscdiff2,guscdiff3,&
+ theta_diff
+
+
+!
+! FP - 30/10/2014 Temporary specifications for homology restraints
+!
+ real(kind=8) :: utheta_i,gutheta_i,sum_gtheta,sum_sgtheta,&
+ sgtheta
+ real(kind=8), dimension (nres) :: guscdiff,usc_diff
+ real(kind=8) :: sum_godl,sgodl,grad_odl3,ggodl,sum_gdih,&
+ sum_guscdiff,sum_sgdih,sgdih,grad_dih3,usc_diff_i,dxx,dyy,dzz,&
+ betai,sum_sgodl,dij,max_template
+! real(kind=8) :: dist,pinorm
+!
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.CONTROL'
+! include 'COMMON.HOMOLOGY'
+! include 'COMMON.QRESTR'
+!
+! From subroutine Econstr_back
+!
+! include 'COMMON.NAMES'
+! include 'COMMON.TIME1'
+!
+
+
+ do i=1,max_template
+ distancek(i)=9999999.9
+ enddo
+
+
+ odleg=0.0d0
+
+! Pseudo-energy and gradient from homology restraints (MODELLER-like
+! function)
+! AL 5/2/14 - Introduce list of restraints
+! write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs start -------"
+#endif
+ do ii = link_start_homo,link_end_homo
+ i = ires_homo(ii)
+ j = jres_homo(ii)
+ dij=dist(i,j)
+! write (iout,*) "dij(",i,j,") =",dij
+ nexl=0
+ do k=1,constr_homology
+! write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
+ if(.not.l_homo(k,ii)) then
+ nexl=nexl+1
+ cycle
+ endif
+ distance(k)=odl(k,ii)-dij
+! write (iout,*) "distance(",k,") =",distance(k)
+!
+! For Gaussian-type Urestr
+!
+ distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument
+! write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii)
+! write (iout,*) "distancek(",k,") =",distancek(k)
+! distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
+!
+! For Lorentzian-type Urestr
+!
+ if (waga_dist.lt.0.0d0) then
+ sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii))
+ distancek(k)=distance(k)**2/(sigma_odlir(k,ii)* &
+ (distance(k)**2+sigma_odlir(k,ii)**2))
+ endif
+ enddo
+
+! min_odl=minval(distancek)
+ if (nexl.gt.0) then
+ min_odl=0.0d0
+ else
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if (l_homo(kk,ii) .and. distancek(kk).lt.min_odl) &
+ min_odl=distancek(kk)
+ enddo
+ endif
+
+! write (iout,* )"min_odl",min_odl
+#ifdef DEBUG
+ write (iout,*) "ij dij",i,j,dij
+ write (iout,*) "distance",(distance(k),k=1,constr_homology)
+ write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
+ write (iout,* )"min_odl",min_odl
+#endif
+#ifdef OLDRESTR
+ odleg2=0.0d0
+#else
+ if (waga_dist.ge.0.0d0) then
+ odleg2=nexl
+ else
+ odleg2=0.0d0
+ endif
+#endif
+ do k=1,constr_homology
+! Nie wiem po co to liczycie jeszcze raz!
+! odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
+! & (2*(sigma_odl(i,j,k))**2))
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ godl(k)=dexp(-distancek(k)+min_odl)
+ odleg2=odleg2+godl(k)
+!
+! For Lorentzian-type Urestr
+!
+ else
+ odleg2=odleg2+distancek(k)
+ endif
+
+!cc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
+!cc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
+!cc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1),
+!cc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
+
+ enddo
+! write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+! write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#endif
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+!
+! For Lorentzian-type Urestr
+!
+ else
+ odleg=odleg+odleg2/constr_homology
+ endif
+!
+! write (iout,*) "odleg",odleg ! sum of -ln-s
+! Gradient
+!
+! For Gaussian-type Urestr
+!
+ if (waga_dist.ge.0.0d0) sum_godl=odleg2
+ sum_sgodl=0.0d0
+ do k=1,constr_homology
+! godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+! & *waga_dist)+min_odl
+! sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+!
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+! For Gaussian-type Urestr
+!
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+!
+! For Lorentzian-type Urestr
+!
+ else
+ sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+ &
+ sigma_odlir(k,ii)**2)**2)
+ endif
+ sum_sgodl=sum_sgodl+sgodl
+
+! sgodl2=sgodl2+sgodl
+! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+! write(iout,*) "constr_homology=",constr_homology
+! write(iout,*) i, j, k, "TEST K"
+ enddo
+! print *, "ok",iset
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ grad_odl3=waga_homology(iset)*waga_dist &
+ *sum_sgodl/(sum_godl*dij)
+! print *, "ok"
+!
+! For Lorentzian-type Urestr
+!
+ else
+! Original grad expr modified by analogy w Gaussian-type Urestr grad
+! grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl
+ grad_odl3=-waga_homology(iset)*waga_dist* &
+ sum_sgodl/(constr_homology*dij)
+! print *, "ok2"
+ endif
+!
+! grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+! write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+! & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+!cc write(iout,*) godl, sgodl, grad_odl3
+
+! grad_odl=grad_odl+grad_odl3
+
+ do jik=1,3
+ ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+!cc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+!cc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl,
+!cc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+ ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+ ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+!cc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+!cc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+! if (i.eq.25.and.j.eq.27) then
+! write(iout,*) "jik",jik,"i",i,"j",j
+! write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl
+! write(iout,*) "grad_odl3",grad_odl3
+! write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j)
+! write(iout,*) "ggodl",ggodl
+! write(iout,*) "ghpbc(",jik,i,")",
+! & ghpbc(jik,i),"ghpbc(",jik,j,")",
+! & ghpbc(jik,j)
+! endif
+ enddo
+!cc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
+!cc & dLOG(odleg2),"-odleg=", -odleg
+
+ enddo ! ii-loop for dist
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs end -------"
+! if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or.
+! & waga_d.eq.1.0d0) call sum_gradient
+#endif
+! Pseudo-energy and gradient from dihedral-angle restraints from
+! homology templates
+! write (iout,*) "End of distance loop"
+! call flush(iout)
+ kat=0.0d0
+! write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs start -------"
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg)
+ enddo
+#endif
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ kat2=0.0d0
+! betai=beta(i,i+1,i+2,i+3)
+ betai = phi(i)
+! write (iout,*) "betai =",betai
+ do k=1,constr_homology
+ dih_diff(k)=pinorm(dih(k,i)-betai)
+!d write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k)
+!d & ,sigma_dih(k,i)
+! if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
+! & -(6.28318-dih_diff(i,k))
+! if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
+! & 6.28318+dih_diff(i,k)
+#ifdef OLD_DIHED
+ kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+ kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#endif
+! kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
+ gdih(k)=dexp(kat3)
+ kat2=kat2+gdih(k)
+! write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+! write(*,*)""
+ enddo
+! write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
+! write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "i",i," betai",betai," kat2",kat2
+ write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+ if (kat2.le.1.0d-14) cycle
+ kat=kat-dLOG(kat2/constr_homology)
+! write (iout,*) "kat",kat ! sum of -ln-s
+
+!cc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+!cc & dLOG(kat2), "-kat=", -kat
+
+! ----------------------------------------------------------------------
+! Gradient
+! ----------------------------------------------------------------------
+
+ sum_gdih=kat2
+ sum_sgdih=0.0d0
+ do k=1,constr_homology
+#ifdef OLD_DIHED
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+#else
+ sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd
+#endif
+! sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
+ sum_sgdih=sum_sgdih+sgdih
+ enddo
+! grad_dih3=sum_sgdih/sum_gdih
+ grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih
+! print *, "ok3"
+
+! write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
+!cc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
+!cc & gloc(nphi+i-3,icg)
+ gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
+! if (i.eq.25) then
+! write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
+! endif
+!cc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+!cc & gloc(nphi+i-3,icg)
+
+ enddo ! i-loop for dih
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs end -------"
+#endif
+
+! Pseudo-energy and gradient for theta angle restraints from
+! homology templates
+! FP 01/15 - inserted from econstr_local_test.F, loop structure
+! adapted
+
+!
+! For constr_homology reference structures (FP)
+!
+! Uconst_back_tot=0.0d0
+ Eval=0.0d0
+ Erot=0.0d0
+! Econstr_back legacy
+ do i=1,nres
+! do i=ithet_start,ithet_end
+ dutheta(i)=0.0d0
+ enddo
+! do i=loc_start,loc_end
+ do i=-1,nres
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+!
+! do iref=1,nref
+! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+! write (iout,*) "waga_theta",waga_theta
+ if (waga_theta.gt.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "usampl",usampl
+ write(iout,*) "------- theta restrs start -------"
+! do i=ithet_start,ithet_end
+! write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg)
+! enddo
+#endif
+! write (iout,*) "maxres",maxres,"nres",nres
+
+ do i=ithet_start,ithet_end
+!
+! do i=1,nfrag_back
+! ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+!
+! Deviation of theta angles wrt constr_homology ref structures
+!
+ utheta_i=0.0d0 ! argument of Gaussian for single k
+ gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+! do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
+! over residues in a fragment
+! write (iout,*) "theta(",i,")=",theta(i)
+ do k=1,constr_homology
+!
+! dtheta_i=theta(j)-thetaref(j,iref)
+! dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
+ theta_diff(k)=thetatpl(k,i)-theta(i)
+!d write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k)
+!d & ,sigma_theta(k,i)
+
+!
+ utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
+! utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
+ gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
+ gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk)
+! Gradient for single Gaussian restraint in subr Econstr_back
+! dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
+!
+ enddo
+! write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps
+! write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps
+
+!
+! Gradient for multiple Gaussian restraint
+ sum_gtheta=gutheta_i
+ sum_sgtheta=0.0d0
+ do k=1,constr_homology
+! New generalized expr for multiple Gaussian from Econstr_back
+ sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd
+!
+! sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
+ sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
+ enddo
+! Final value of gradient using same var as in Econstr_back
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) &
+ +sum_sgtheta/sum_gtheta*waga_theta &
+ *waga_homology(iset)
+! print *, "ok4"
+
+! dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+! & *waga_homology(iset)
+! dutheta(i)=sum_sgtheta/sum_gtheta
+!
+! Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
+ Eval=Eval-dLOG(gutheta_i/constr_homology)
+! write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps
+! write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s
+! Uconst_back=Uconst_back+utheta(i)
+ enddo ! (i-loop for theta)
+#ifdef DEBUG
+ write(iout,*) "------- theta restrs end -------"
+#endif
+ endif
+!
+! Deviation of local SC geometry
+!
+! Separation of two i-loops (instructed by AL - 11/3/2014)
+!
+! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+! write (iout,*) "waga_d",waga_d
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs start -------"
+ write (iout,*) "Initial duscdiff,duscdiffx"
+ do i=loc_start,loc_end
+ write (iout,*) i,(duscdiff(jik,i),jik=1,3), &
+ (duscdiffx(jik,i),jik=1,3)
+ enddo
+#endif
+ do i=loc_start,loc_end
+ usc_diff_i=0.0d0 ! argument of Gaussian for single k
+ guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+! do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+! write(iout,*) "xxtab, yytab, zztab"
+! write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+!
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+! Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+! write(iout,*) "dxx, dyy, dzz"
+!d write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i)
+!
+ usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
+! usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
+! uscdiffk(k)=usc_diff(i)
+ guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
+! write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i),
+! & " guscdiff2",guscdiff2(k)
+ guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk)
+! write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+! & xxref(j),yyref(j),zzref(j)
+ enddo
+!
+! Gradient
+!
+! Generalized expression for multiple Gaussian acc to that for a single
+! Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+!
+! Original implementation
+! sum_guscdiff=guscdiff(i)
+!
+! sum_sguscdiff=0.0d0
+! do k=1,constr_homology
+! sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+! sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+! sum_sguscdiff=sum_sguscdiff+sguscdiff
+! enddo
+!
+! Implementation of new expressions for gradient (Jan. 2015)
+!
+! grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+ do k=1,constr_homology
+!
+! New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+! before. Now the drivatives should be correct
+!
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+! Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+ sum_guscdiff=guscdiff2(k)* &!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ sigma_d(k,i) ! for the grad wrt r'
+! sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+
+!
+! New implementation
+ sum_guscdiff = waga_homology(iset)*waga_d*sum_guscdiff
+ do jik=1,3
+ duscdiff(jik,i-1)=duscdiff(jik,i-1)+ &
+ sum_guscdiff*(dXX_C1tab(jik,i)*dxx+ &
+ dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i)
+ duscdiff(jik,i)=duscdiff(jik,i)+ &
+ sum_guscdiff*(dXX_Ctab(jik,i)*dxx+ &
+ dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i)
+ duscdiffx(jik,i)=duscdiffx(jik,i)+ &
+ sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+ &
+ dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i)
+! print *, "ok5"
+!
+#ifdef DEBUG
+! write(iout,*) "jik",jik,"i",i
+ write(iout,*) "dxx, dyy, dzz"
+ write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+ write(iout,*) "guscdiff2(",k,")",guscdiff2(k)
+ write(iout,*) "sum_sguscdiff",sum_guscdiff,waga_homology(iset),waga_d
+ write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+ write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+ write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+ write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+ write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+ write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+ write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+ write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+ write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+ write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+ write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+ write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+! endif
+#endif
+ enddo
+ enddo
+! print *, "ok6"
+!
+! uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+! usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+!
+! write (iout,*) i," uscdiff",uscdiff(i)
+!
+! Put together deviations from local geometry
+
+! Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
+! & wfrag_back(3,i,iset)*uscdiff(i)
+ Erot=Erot-dLOG(guscdiff(i)/constr_homology)
+! write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps
+! write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s
+! Uconst_back=Uconst_back+usc_diff(i)
+!
+! Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?)
+!
+! New implment: multiplied by sum_sguscdiff
+!
+
+ enddo ! (i-loop for dscdiff)
+
+! endif
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs end -------"
+ write (iout,*) "------ After SC loop in e_modeller ------"
+ do i=loc_start,loc_end
+ write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+ write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
+ enddo
+ if (waga_theta.eq.1.0d0) then
+ write (iout,*) "in e_modeller after SC restr end: dutheta"
+ do i=ithet_start,ithet_end
+ write (iout,*) i,dutheta(i)
+ enddo
+ endif
+ if (waga_d.eq.1.0d0) then
+ write (iout,*) "e_modeller after SC loop: duscdiff/x"
+ do i=1,nres
+ write (iout,*) i,(duscdiff(j,i),j=1,3)
+ write (iout,*) i,(duscdiffx(j,i),j=1,3)
+ enddo
+ endif
+#endif
+
+! Total energy from homology restraints
+#ifdef DEBUG
+ write (iout,*) "odleg",odleg," kat",kat
+#endif
+!
+! Addition of energy of theta angle and SC local geom over constr_homologs ref strs
+!
+! ehomology_constr=odleg+kat
+!
+! For Lorentzian-type Urestr
+!
+
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ ehomology_constr=(waga_dist*odleg+waga_angle*kat+ &
+ waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+! write (iout,*) "ehomology_constr=",ehomology_constr
+! print *, "ok7"
+ else
+!
+! For Lorentzian-type Urestr
+!
+ ehomology_constr=(-waga_dist*odleg+waga_angle*kat+ &
+ waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+! write (iout,*) "ehomology_constr=",ehomology_constr
+ print *, "ok8"
+ endif
+#ifdef DEBUG
+ write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, &
+ "Eval",waga_theta,eval, &
+ "Erot",waga_d,Erot
+ write (iout,*) "ehomology_constr",ehomology_constr
+#endif
+ return
+!
+! FP 01/15 end
+!
+ 748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+ 747 format(a12,i4,i4,i4,f8.3,f8.3)
+ 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+ 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+ 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X, &
+ f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
+ end subroutine e_modeller
+
+!----------------------------------------------------------------------------
subroutine ebend_kcc(etheta)
logical lprn
double precision thybt1(maxang_kcc),etheta
wscbase*gvdwc_scbase(j,i)+ &
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+ &
- wpeppho*gvdwc_peppho(j,i)
+ wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)
wscbase*gvdwc_scbase(j,i)+ &
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+&
- wpeppho*gvdwc_peppho(j,i)
+ wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)
enddo
+wbond_nucl*gradb_nucl(j,i) &
+0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)&
+wvdwpsb*gvdwpsb1(j,i))&
- +wsbloc*gsbloc(j,i)
+ +wsbloc*gsbloc(j,i)+wcatnucl*gradnuclcat(j,i)
+wcatprot* gradpepcatx(j,i)&
+wscbase*gvdwx_scbase(j,i) &
+wpepbase*gvdwx_pepbase(j,i)&
- +wscpho*gvdwx_scpho(j,i)
+ +wscpho*gvdwx_scpho(j,i)+wcatnucl*gradnuclcatx(j,i)
! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
enddo
enddo
+! write(iout,*), "const_homol",constr_homology
+ if (constr_homology.gt.0) then
+ do i=1,nct
+ do j=1,3
+ gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
+! write(iout,*) "duscdiff",duscdiff(j,i)
+ gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i)
+ enddo
+ enddo
+ endif
!#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc before adding corr"
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)
subroutine check_ecartint
! Check the gradient of the energy in Cartesian coordinates.
use io_base, only: intout
+ use MD_data, only: iset
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.CONTROL'
icg=1
nf=0
nfl=0
+ if (iset.eq.0) iset=1
call intout
! call intcartderiv
! call checkintcartgrad
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
subroutine check_ecartint
! Check the gradient of the energy in Cartesian coordinates.
use io_base, only: intout
+ use MD_data, only: iset
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.CONTROL'
icg=1
nf=0
nfl=0
+ if (iset.eq.0) iset=1
call intout
! call intcartderiv
! call checkintcartgrad
enddo
do j=1,3
grad_s(j,0)=gcart(j,0)
+ grad_s(j+3,0)=gxcart(j,0)
enddo
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
integer :: i,n_corr,n_corr1,ierror,ierr
real(kind=8) :: evdw2,evdw2_14,ehpb,etors,edihcnstr,etors_d,esccor,&
evdw,ees,evdw1,eel_loc,eello_turn3,eello_turn4,&
- ecorr,ecorr5,ecorr6,eturn6,time00
+ ecorr,ecorr5,ecorr6,eturn6,time00, ehomology_constr
! write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
!elwrite(iout,*)"in etotal long"
#endif
endif
!elwrite(iout,*)"in etotal long"
-
+ ehomology_constr=0.0d0
#ifdef MPI
! write(iout,*) "ETOTAL_LONG Processor",fg_rank,
! & " absolute rank",myrank," nfgtasks",nfgtasks
energia(9)=eello_turn4
energia(10)=eturn6
energia(20)=Uconst+Uconst_back
+ energia(51)=ehomology_constr
call sum_energy(energia,.true.)
! write (iout,*) "Exit ETOTAL_LONG"
call flush(iout)
!el local variables
integer :: i,nres6
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors
- real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr
+ real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr, &
+ ehomology_constr
nres6=6*nres
! write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
call etor_d(etors_d)
endif
!
+! Homology restraints
+!
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+! print *,"tu"
+ else
+ ehomology_constr=0.0d0
+ endif
+
+!
! 21/5/07 Calculate local sicdechain correlation energy
!
if (wsccor.gt.0.0d0) then
energia(17)=estr
energia(19)=edihcnstr
energia(21)=esccor
+ energia(51)=ehomology_constr
! write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
call flush(iout)
call sum_energy(energia,.true.)
#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
gvdwx_scpho(j,i)=0.0d0
gvdwc_scpho(j,i)=0.0d0
gvdwc_peppho(j,i)=0.0d0
+ gradnuclcatx(j,i)=0.0d0
+ gradnuclcat(j,i)=0.0d0
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
enddo
enddo
do i=0,nres
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
else
maxconts=10*nres ! (maxconts=maxres/4)
endif
- maxcont=12*nres ! Max. number of SC contacts
+ maxcont=100*nres ! Max. number of SC contacts
maxvar=6*nres ! Max. number of variables
!el maxdim=(nres-1)*(nres-2)/2 ! Max. number of derivatives of virtual-bond
maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond
allocate(gacontm_hb3(3,maxconts,nres))
allocate(gacont_hbr(3,maxconts,nres))
allocate(grij_hb_cont(3,maxconts,nres))
-!(3,maxconts,maxres)
+ !(3,maxconts,maxres)
allocate(facont_hb(maxconts,nres))
allocate(ees0p(maxconts,nres))
allocate(gradpepcat(3,-1:nres))
allocate(gradpepcatx(3,-1:nres))
allocate(gradcatcat(3,-1:nres))
+ allocate(gradnuclcat(3,-1:nres))
+ allocate(gradnuclcatx(3,-1:nres))
!(3,maxres)
allocate(grad_shield_side(3,maxcontsshi,-1:nres))
allocate(grad_shield_loc(3,maxcontsshi,-1:nres))
allocate(dutheta(nres))
allocate(dugamma(nres))
!(maxres)
- allocate(duscdiff(3,nres))
- allocate(duscdiffx(3,nres))
+ allocate(duscdiff(3,-1:nres))
+ allocate(duscdiffx(3,-1:nres))
!(3,maxres)
!el i io:read_fragments
! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)
!(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
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(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
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
- +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
yj=c(2,j)
zj=c(3,j)
call to_box(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
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
- +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
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)
- 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
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
- +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_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
+! write(iout,*) "KURWA0",d1
CALL edq_cat(ecl, elj, epol)
eheadtail = 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
end subroutine ecat_prot
!----------------------------------------------------------------------------
+!---------------------------------------------------------------------------
+ subroutine ecat_nucl(ecation_nucl)
+ integer i,j,k,subchap,itmp,inum,itypi,itypj
+ real(kind=8) :: xi,yi,zi,xj,yj,zj
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_nucl,Evan1,Evan2,Ecav,Egb,wdip1,wdip2, &
+ wvan1,wvan2,wgbsig,wgbeps,wgbchi,wgbchip,wcav1,wcav2,wcav3,wcav4, &
+ wcavsig,wcavchi,v1m,v1dpdx,wh2o,wc,Edip,rcs2,invrcs6,invrcs8,invrcs12, &
+ invrcs14,rcb,rcb2,invrcb,invrcb2,invrcb4,invrcb6,cosinus,cos2,dcosdcatconst, &
+ dcosdcalpconst,dcosdcmconst,rcav,rcav11,rcav12,constcav1,constcav2, &
+ constgb1,constgb2,constdvan1,constdvan2,sgb,sgb6,sgb7,sgb12,sgb13, &
+ cavnum,cavdenom,invcavdenom2,dcavnumdcos,dcavnumdr,dcavdenomdcos, &
+ dcavdenomdr,sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+ real(kind=8),dimension(3) ::gg,r,dEtotalCm,dEtotalCalp,dEvan1Cm,&
+ 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,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
+ itmp=itmp+nres_molec(i)
+ enddo
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ if ((itype(i,2).eq.ntyp1_molec(2))) cycle ! leave dummy atoms
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ do j=itmp+1,itmp+nres_molec(5)
+ xj=c(1,j)
+ 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_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)
+ itypj=itype(j,5)
+ do k=1,13
+ vcatnuclprm(k)=catnuclprm(k,itypi,itypj)
+ enddo
+ do k=1,3
+ vcm(k)=c(k,i+nres)
+ 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
+ dx(k)=boxshift(vcat(k)-vcm(k),boxik(k))
+! do k=1,3
+ v1(k)=dc(k,i+nres)
+ 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)
+! The weights of the energy function calculated from
+!The quantum mechanical Gaussian simulations of potassium and sodium with deoxynucleosides
+ wh2o=78
+ wdip1 = vcatnuclprm(1)
+ wdip1 = wdip1/wh2o !w1
+ wdip2 = vcatnuclprm(2)
+ wdip2 = wdip2/wh2o !w2
+ wvan1 = vcatnuclprm(3)
+ wvan2 = vcatnuclprm(4) !pis1
+ wgbsig = vcatnuclprm(5) !sigma0
+ wgbeps = vcatnuclprm(6) !epsi0
+ wgbchi = vcatnuclprm(7) !chi1
+ wgbchip = vcatnuclprm(8) !chip1
+ wcavsig = vcatnuclprm(9) !sig
+ wcav1 = vcatnuclprm(10) !b1
+ wcav2 = vcatnuclprm(11) !b2
+ wcav3 = vcatnuclprm(12) !b3
+ wcav4 = vcatnuclprm(13) !b4
+ wcavchi = vcatnuclprm(14) !chis1
+ rcs2 = v2(1)**2+v2(2)**2+v2(3)**2
+ invrcs6 = 1/rcs2**3
+ invrcs8 = invrcs6/rcs2
+ invrcs12 = invrcs6**2
+ invrcs14 = invrcs12/rcs2
+ rcb2 = dx(1)**2+dx(2)**2+dx(3)**2
+ rcb = sqrt(rcb2)
+ invrcb = 1/rcb
+ invrcb2 = invrcb**2
+ invrcb4 = invrcb2**2
+ invrcb6 = invrcb4*invrcb2
+ cosinus = v1dpdx/(v1m*rcb)
+ cos2 = cosinus**2
+ dcosdcatconst = invrcb2/v1m
+ dcosdcalpconst = invrcb/v1m**2
+ dcosdcmconst = invrcb2/v1m**2
+ do k=1,3
+ dcosdcat(k) = (v1(k)*rcb-dx(k)*v1m*cosinus)*dcosdcatconst
+ dcosdcalp(k) = (v1(k)*rcb*cosinus-dx(k)*v1m)*dcosdcalpconst
+ dcosdcm(k) = ((dx(k)-v1(k))*v1m*rcb+ &
+ cosinus*(dx(k)*v1m**2-v1(k)*rcb2))*dcosdcmconst
+ enddo
+ rcav = rcb/wcavsig
+ rcav11 = rcav**11
+ rcav12 = rcav11*rcav
+ constcav1 = 1-wcavchi*cos2
+ constcav2 = sqrt(constcav1)
+ constgb1 = 1/sqrt(1-wgbchi*cos2)
+ constgb2 = wgbeps*(1-wgbchip*cos2)**2
+ constdvan1 = 12*wvan1*wvan2**12*invrcs14
+ constdvan2 = 6*wvan1*wvan2**6*invrcs8
+!----------------------------------------------------------------------------
+!Gay-Berne term
+!---------------------------------------------------------------------------
+ sgb = 1/(1-constgb1+(rcb/wgbsig))
+ sgb6 = sgb**6
+ sgb7 = sgb6*sgb
+ sgb12 = sgb6**2
+ sgb13 = sgb12*sgb
+ Egb = constgb2*(sgb12-sgb6)
+ do k=1,3
+ dEgbdCat(k) = -constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcat(k)
+ dEgbdCm(k) = constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcm(k)
+ dEgbdCalp(k) = (constgb1**3*constgb2*wgbchi*cosinus &
+ *(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcalp(k)
+ enddo
+!----------------------------------------------------------------------------
+!cavity term
+!---------------------------------------------------------------------------
+ cavnum = sqrt(rcav*constcav2)+wcav2*rcav*constcav2-wcav3
+ cavdenom = 1+wcav4*rcav12*constcav1**6
+ Ecav = wcav1*cavnum/cavdenom
+ invcavdenom2 = 1/cavdenom**2
+ dcavnumdcos = -wcavchi*cosinus/constcav2 &
+ *(sqrt(rcav/constcav2)/2+wcav2*rcav)
+ dcavnumdr = (0.5*sqrt(constcav2/rcav)+wcav2*constcav2)/wcavsig
+ dcavdenomdcos = -12*wcav4*wcavchi*rcav12*constcav1**5*cosinus
+ dcavdenomdr = 12*wcav4/wcavsig*rcav11*constcav1**6
+ do k=1,3
+ dEcavdCat(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcat(k)+(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCm(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcm(k)-(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCalp(k) = (dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcalp(k)*wcav1*invcavdenom2
+ enddo
+!----------------------------------------------------------------------------
+!van der Waals and dipole-charge interaction energy
+!---------------------------------------------------------------------------
+ Evan1 = wvan1*wvan2**12*invrcs12
+ do k=1,3
+ dEvan1Cat(k) = -v2(k)*constdvan1
+ dEvan1Cm(k) = 0.0d0
+ dEvan1Calp(k) = v2(k)*constdvan1
+ enddo
+ Evan2 = -wvan1*wvan2**6*invrcs6
+ do k=1,3
+ dEvan2Cat(k) = v2(k)*constdvan2
+ dEvan2Cm(k) = 0.0d0
+ dEvan2Calp(k) = -v2(k)*constdvan2
+ enddo
+ Edip = wdip1*cosinus*invrcb2-wdip2*(1-cos2)*invrcb4
+ do k=1,3
+ dEdipCat(k) = (-2*wdip1*cosinus*invrcb4 &
+ +4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcat(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCm(k) = (2*wdip1*cosinus*invrcb4 &
+ -4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcm(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCalp(k) = dcosdcalp(k)*(wdip1*invrcb2 &
+ +2*wdip2*cosinus*invrcb4)
+ enddo
+ if (energy_dec) write (iout,'(2i5,4(a6,f7.3))') i,j, &
+ ' E GB ',Egb,' ECav ',Ecav,' Evdw ',Evan1+Evan2,' Edip ',Edip
+ ecation_nucl=ecation_nucl+Ecav+Egb+Edip+Evan1+Evan2
+ do k=1,3
+ dEtotalCat(k) = dEcavdCat(k)+dEvan1Cat(k)+dEvan2Cat(k) &
+ +dEgbdCat(k)+dEdipCat(k)
+ dEtotalCm(k) = dEcavdCm(k)+dEvan1Cm(k)+dEvan2Cm(k) &
+ +dEgbdCm(k)+dEdipCm(k)
+ dEtotalCalp(k) = dEcavdCalp(k)+dEgbdCalp(k)+dEvan1Calp(k) &
+ +dEdipCalp(k)+dEvan2Calp(k)
+ enddo
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradnuclcatx(k,i)=gradnuclcatx(k,i)+dEtotalCm(k)
+ gradnuclcat(k,i)=gradnuclcat(k,i)+gg(k)
+ gradnuclcat(k,j)=gradnuclcat(k,j)+dEtotalCat(k)
+ enddo
+ enddo !j
+ enddo !i
+ return
+ end subroutine ecat_nucl
+
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
subroutine eprot_sc_base(escbase)
yj=c(2,j+nres)
zj=c(3,j+nres)
call to_box(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
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
- +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
yj=(c(2,j)+c(2,j+1))/2.0
zj=(c(3,j)+c(3,j+1))/2.0
call to_box(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
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
- +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! 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
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
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
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)
!c!-------------------------------------------------------------------
!c! ecl
- write(iout,*) "KURWA2",Rhead
+! write(iout,*) "KURWA2",Rhead
sparrow = w1 * Qj * om1
hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
ECL = sparrow / Rhead**2.0d0 &
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)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dist_init=xj**2+yj**2+zj**2
+! dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
! 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
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
! 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
- write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
+! write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
ilist_pp=0
r_buff_list=5.0
do i=iatel_s,iatel_e
projects / unres4.git / blobdiff