X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=fce83dda9a05671b28767dc22ac3b93a171d262e;hb=41585054482bea08de6120a0ec991bc5288c9fe5;hp=3870bb0aaeef6d24551f54ef0c7b9693d9426464;hpb=da0ffdc0802840f3e36a2f41166ad72ced8f7845;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 3870bb0..fce83dd 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -1,4 +1,4 @@ - module energy + module energy !----------------------------------------------------------------------------- use io_units use names @@ -236,7 +236,7 @@ ! 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, & @@ -247,7 +247,7 @@ ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& ecorr3_nucl ! energies for ions - real(kind=8) :: ecation_prot,ecationcation + real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber ! energies for protein nucleic acid interaction real(kind=8) :: escbase,epepbase,escpho,epeppho @@ -264,8 +264,9 @@ 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 @@ -329,8 +330,9 @@ weights_(41)=wcatcat weights_(42)=wcatprot weights_(46)=wscbase - weights_(47)=wscpho - weights_(48)=wpeppho + weights_(47)=wpepbase + weights_(48)=wscpho + weights_(49)=wpeppho ! wcatcat= weights(41) ! wcatprot=weights(42) @@ -376,13 +378,29 @@ wcatcat= weights(41) wcatprot=weights(42) wscbase=weights(46) - wscpho=weights(47) - wpeppho=weights(48) + wpepbase=weights(47) + wscpho=weights(48) + wpeppho=weights(49) +! welpsb=weights(28)*fact(1) +! +! wcorr_nucl= weights(37)*fact(1) +! wcorr3_nucl=weights(38)*fact(2) +! wtor_nucl= weights(35)*fact(1) +! wtor_d_nucl=weights(36)*fact(2) + endif time_Bcast=time_Bcast+MPI_Wtime()-time00 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 @@ -793,6 +811,8 @@ call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) + else + Eafmforce=0.0d0 endif endif if (tubemode.eq.1) then @@ -821,6 +841,7 @@ 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 @@ -829,9 +850,11 @@ 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" +! 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) @@ -839,8 +862,16 @@ else call ecatcat(ecationcation) endif + if (oldion.gt.0) then call ecat_prot(ecation_prot) - 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) @@ -852,7 +883,7 @@ epeppho=0.0 endif ! call ecatcat(ecationcation) -! print *,"after ebend", ebe_nucl +! print *,"after ebend", wtor_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -916,12 +947,13 @@ ! Here are the energies showed per procesor if the are more processors ! per molecule then we sum it up in sum_energy subroutine ! print *," Processor",myrank," calls SUM_ENERGY" - energia(41)=ecation_prot - energia(42)=ecationcation + energia(42)=ecation_prot + energia(41)=ecationcation energia(46)=escbase energia(47)=epepbase energia(48)=escpho energia(49)=epeppho +! energia(50)=ecations_prot_amber call sum_energy(energia,.true.) if (dyn_ss) call dyn_set_nss ! print *," Processor",myrank," left SUM_ENERGY" @@ -964,7 +996,7 @@ 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 + real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber real(kind=8) :: escbase,epepbase,escpho,epeppho integer :: i #ifdef MPI @@ -1042,12 +1074,14 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) - ecation_prot=energia(41) - ecationcation=energia(42) + ecation_prot=energia(42) + ecationcation=energia(41) escbase=energia(46) epepbase=energia(47) escpho=energia(48) epeppho=energia(49) +! ecations_prot_amber=energia(50) + ! energia(41)=ecation_prot ! energia(42)=ecationcation @@ -1187,7 +1221,12 @@ wtor=weights(13)*fact(1) wtor_d=weights(14)*fact(2) wsccor=weights(21)*fact(1) - + welpsb=weights(28)*fact(1) + wcorr_nucl= weights(37)*fact(1) + wcorr3_nucl=weights(38)*fact(2) + wtor_nucl= weights(35)*fact(1) + wtor_d_nucl=weights(36)*fact(2) + wpepbase=weights(47)*fact(1) return end subroutine rescale_weights !----------------------------------------------------------------------------- @@ -1207,7 +1246,7 @@ 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 + real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber real(kind=8) :: escbase,epepbase,escpho,epeppho etot=energia(0) @@ -1255,12 +1294,13 @@ etors_d_nucl=energia(36) ecorr_nucl=energia(37) ecorr3_nucl=energia(38) - ecation_prot=energia(41) - ecationcation=energia(42) + ecation_prot=energia(42) + ecationcation=energia(41) escbase=energia(46) epepbase=energia(47) escpho=energia(48) epeppho=energia(49) +! ecations_prot_amber=energia(50) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -1331,11 +1371,11 @@ ecorr,wcorr,& 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,Eafmforc, & + ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce, & etube,wtube, & 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& + 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,& @@ -1799,7 +1839,7 @@ ! 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,& @@ -1822,7 +1862,11 @@ 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 @@ -1871,8 +1915,8 @@ ! ! 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 @@ -2013,8 +2057,8 @@ ! write(iout,*)"c ", c(1,:), c(2,:), c(3,:) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) ! print *,sss_ele_cut,sss_ele_grad,& ! 1.0d0/(rij),r_cut_ele,rlamb_ele if (sss_ele_cut.le.0.0) cycle @@ -2076,7 +2120,7 @@ fac=rij*fac ! print *,'before fac',fac,rij,evdwij fac=fac+evdwij*sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij + *rij ! print *,'grad part scale',fac, & ! evdwij*sss_ele_grad/sss_ele_cut & ! /sigma(itypi,itypj)*rij @@ -2097,8 +2141,8 @@ ! 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 @@ -2837,7 +2881,17 @@ #endif #else if (i.gt. nnt+2 .and. i.lt.nct+2) then +! 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 else iti=nloctyp endif @@ -2848,6 +2902,7 @@ else iti1=nloctyp endif +! print *,i,iti b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) @@ -3012,7 +3067,7 @@ ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then if (itype(i-1,1).eq.0) then - iti1=ntortyp+1 + iti1=nloctyp elseif (itype(i-1,1).le.ntyp) then iti1 = itype2loc(itype(i-1,1)) else @@ -3037,7 +3092,7 @@ call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2)) call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2)) ! Vectors and matrices dependent on a single virtual-bond dihedral. - call matvec2(DD(1,1,i-2),b1tilde(1,iti1),auxvec(1)) + call matvec2(DD(1,1,i-2),b1tilde(1,i-1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2)) @@ -3376,7 +3431,7 @@ 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 @@ -3576,7 +3631,11 @@ ! 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) @@ -3618,11 +3677,11 @@ ! 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 @@ -3796,7 +3855,7 @@ ! 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 @@ -4012,11 +4071,11 @@ !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 @@ -4345,7 +4404,9 @@ +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), @@ -4363,7 +4424,8 @@ 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 @@ -4378,7 +4440,8 @@ 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=& @@ -4392,7 +4455,9 @@ 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 @@ -4591,10 +4656,12 @@ 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 @@ -4673,28 +4740,36 @@ 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! @@ -4841,10 +4916,15 @@ !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 @@ -5053,9 +5133,9 @@ a_temp(1,2)=a23 a_temp(2,1)=a32 a_temp(2,2)=a33 - iti1=itortyp(itype(i+1,1)) - iti2=itortyp(itype(i+2,1)) - iti3=itortyp(itype(i+3,1)) + iti1=i+1 + iti2=i+2 + iti3=i+3 ! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) @@ -5086,7 +5166,7 @@ call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) !c auxilary matrix auxgEvec1 of E matix with Ub2 constant call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) !c derivative of theta i+1 with constant i+3 gs13=scalar2(gtb1(1,i+1),auxvec(1)) !c derivative of theta i+2 with constant i+1 @@ -5186,7 +5266,7 @@ call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) & @@ -5482,7 +5562,7 @@ ! 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,& @@ -5493,7 +5573,10 @@ 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)) @@ -5506,9 +5589,9 @@ 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 @@ -5620,9 +5703,9 @@ 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 @@ -6913,6 +6996,8 @@ ! & 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 @@ -11320,7 +11405,7 @@ +wcorr3_nucl*gradcorr3_nucl(j,i) +& wcatprot* gradpepcat(j,i)+ & wcatcat*gradcatcat(j,i)+ & - wscbase*gvdwc_scbase(j,i) & + wscbase*gvdwc_scbase(j,i)+ & wpepbase*gvdwc_pepbase(j,i)+& wscpho*gvdwc_scpho(j,i)+& wpeppho*gvdwc_peppho(j,i) @@ -11551,7 +11636,7 @@ +gradafm(j,i) & +wliptran*gliptranc(j,i) & +welec*gshieldc(j,i) & - +welec*gshieldc_loc(j,) & + +welec*gshieldc_loc(j,i) & +wcorr*gshieldc_ec(j,i) & +wcorr*gshieldc_loc_ec(j,i) & +wturn3*gshieldc_t3(j,i) & @@ -11640,7 +11725,7 @@ 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 @@ -11656,14 +11741,14 @@ 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 @@ -11863,6 +11948,90 @@ enddo return end subroutine sc_grad + + subroutine sc_grad_cat + 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,nres+i)-om1*erij(k)) + dcosom2(k)=rij*(dc_norm(k,j)-om2*erij(k)) + enddo + do k=1,3 + gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)) +!C print *,'gg',k,gg(k) + enddo +! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut +! write (iout,*) "gg",(gg(k),k=1,3) + do k=1,3 + 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 + +! 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 +! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & +! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + enddo +! +! Calculate the components of the gradient in DC and X +! + do l=1,3 + 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 !----------------------------------------------------------------------------- subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t) @@ -12536,7 +12705,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-4 + aincr=1.0D-5 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 @@ -12763,7 +12932,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-7 + aincr=1.0D-6 write(iout,*) 'Calling CHECK_ECARTINT.',aincr nf=0 icall=0 @@ -14019,8 +14188,8 @@ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj))) if (sss_ele_cut.le.0.0) cycle if (sss.lt.1.0d0) then @@ -14076,7 +14245,7 @@ sigder=fac*sigder fac=rij*fac fac=fac+evdwij*(sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij & + *rij-sss_grad/(1.0-sss)*rij & /sigmaii(itypi,itypj)) ! fac=0.0d0 ! Calculate the radial part of the gradient @@ -14310,8 +14479,8 @@ rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj))) sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj))) - sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj))) - sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj))) + sss_ele_cut=sscale_ele(1.0d0/(rij)) + sss_ele_grad=sscagrad_ele(1.0d0/(rij)) if (sss_ele_cut.le.0.0) cycle if (sss.gt.0.0d0) then @@ -14367,7 +14536,7 @@ sigder=fac*sigder fac=rij*fac fac=fac+evdwij*(sss_ele_grad/sss_ele_cut& - /sigma(itypi,itypj)*rij+sss_grad/sss*rij & + *rij+sss_grad/sss*rij & /sigmaii(itypi,itypj)) ! fac=0.0d0 @@ -17001,7 +17170,7 @@ !#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) @@ -17355,6 +17524,10 @@ 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 @@ -17381,7 +17554,7 @@ #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)) @@ -20396,11 +20569,11 @@ integer :: i,j if(nres.lt.100) then - maxconts=nres + maxconts=10*nres elseif(nres.lt.200) then - maxconts=0.8*nres ! Max. number of contacts per residue + maxconts=10*nres ! Max. number of contacts per residue else - maxconts=0.6*nres ! (maxconts=maxres/4) + maxconts=10*nres ! (maxconts=maxres/4) endif maxcont=12*nres ! Max. number of SC contacts maxvar=6*nres ! Max. number of variables @@ -20828,6 +21001,13 @@ 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 @@ -21775,13 +21955,13 @@ enddo ! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and. IF ( j.gt.i+1 .and.& - num_conti.le.maxconts) THEN + num_conti.le.maxcont) THEN !C !C Calculate the contact function. The ith column of the array JCONT will !C contain the numbers of atoms that make contacts with the atom I (of numbers !C greater than I). The arrays FACONT and GACONT will contain the values of !C the contact function and its derivative. - r0ij=2.20D0*sigma(itypi,itypj) + r0ij=2.20D0*sigma_nucl(itypi,itypj) !c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont) !c write (2,*) "fcont",fcont @@ -21791,7 +21971,7 @@ if (num_conti.gt.maxconts) then write (iout,*) 'WARNING - max. # of contacts exceeded;',& - ' will skip next contacts for this conf.' + ' will skip next contacts for this conf.',maxconts else jcont_hb(num_conti,i)=j !c write (iout,*) "num_conti",num_conti, @@ -22544,7 +22724,7 @@ !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, & @@ -22558,7 +22738,7 @@ 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 @@ -22570,7 +22750,8 @@ 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) @@ -22579,6 +22760,9 @@ 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) @@ -22628,8 +22812,8 @@ ! 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 @@ -22646,7 +22830,8 @@ 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 @@ -22654,96 +22839,812 @@ return end subroutine ecatcat !--------------------------------------------------------------------------- - subroutine ecat_prot(ecation_prot) - 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 +! new for K+ + subroutine ecats_prot_amber(evdw) +! subroutine ecat_prot2(ecation_prot) + use calc_data + use comm_momo + + logical :: lprn +!el local variables + integer :: iint,itypi1,subchap,isel,itmp + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi + real(kind=8) :: evdw + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init,ssgradlipi,ssgradlipj, & + sslipi,sslipj,faclip,alpha_sco + integer :: ii + real(kind=8) :: fracinbuf + real (kind=8) :: escpho + real (kind=8),dimension(4):: ener + real(kind=8) :: b1,b2,egb + real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,& + Lambf,& + Chif,ChiLambf,Fcav,dFdR,dFdOM1,& + ecations_prot_amber,dFdOM2,dFdL,dFdOM12,& + federmaus,& + d1i,d1j +! real(kind=8),dimension(3,2)::erhead_tail +! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance + real(kind=8) :: facd4, adler, Fgb, facd3 + integer troll,jj,istate + real (kind=8) :: dcosom1(3),dcosom2(3) + + evdw=0.0D0 + if (nres_molec(5).eq.0) return + eps_out=80.0d0 +! sss_ele_cut=1.0d0 + itmp=0 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 -! 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=iabs(itype(i,1)) + +! if (i.ne.47) 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) + zi=c(3,nres+i) + 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,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=vbld_inv(i+nres) 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(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=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 + 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 + +! dxj = dc_norm( 1, nres+j ) +! dyj = dc_norm( 2, nres+j ) +! dzj = dc_norm( 3, nres+j ) + + itypi = itype(i,1) + 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+nres) + 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) +! 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+nres) + d1 * dc_norm(k, i+nres) + 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 + 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+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) + 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)) +! gvdwx(k,j) = gvdwx(k,j) & +! + (( dFdR + gg(k) ) * pom) + gradpepcat(k,i) = gradpepcat(k,i) & + - (( dFdR + gg(k) ) * ertail(k)) + 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) + 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) 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 + 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 +! eheadtail = 0.0d0 + + 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 + Qij=Qij*2 + endif + if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then + Qj=Qj*2 + Qij=Qij*2 + endif + write(iout,*) "KURWA0",d1 + + CALL edq_cat(ecl, elj, epol) + eheadtail = ECL + elj + epol +! eheadtail = 0.0d0 + + 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 + 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 +! eheadtail = 0.0d0 + +! ELSE IF ((isel.eq.2.and. & +! iabs(Qi).eq.1).and. & +! nstate(itypi,itypj).ne.1) THEN +!c! Different charge-charge interaction ( +/- or -/+ ) +! 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 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 + + 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 +! 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 + 17 continue + do i=ibond_start,ibond_end + +! 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=dmod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + 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) + +! 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=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 + + 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 @@ -25004,6 +25905,7 @@ 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 @@ -25513,7 +26415,7 @@ use calc_data use comm_momo real (kind=8) :: facd3, facd4, federmaus, adler,& - Ecl,Egb,Epol,Fisocav,Elj,Fgb + Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap ! integer :: k !c! Epol and Gpol analytical parameters alphapol1 = alphapol(itypi,itypj) @@ -25558,10 +26460,15 @@ dGCLdOM12 = 0.0d0 ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq)) Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0) - Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb + debkap=debaykap(itypi,itypj) + Egb = -(332.0d0 * Qij *& + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb ! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out !c! Derivative of Egb is Ggb... - dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb) + dGGBdFGB = -(-332.0d0 * Qij * & + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)& + -(332.0d0 * Qij *& + (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) dGGBdR = dGGBdFGB * dFGBdR !c!------------------------------------------------------------------- @@ -25693,6 +26600,196 @@ END DO RETURN END SUBROUTINE eqq + + SUBROUTINE eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) + use calc_data + use comm_momo + real (kind=8) :: facd3, facd4, federmaus, adler,& + Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap +! integer :: k +!c! Epol and Gpol analytical parameters + alphapol1 = alphapolcat(itypi,itypj) + alphapol2 = alphapolcat(itypj,itypi) +!c! Fisocav and Gisocav analytical parameters + al1 = alphisocat(1,itypi,itypj) + al2 = alphisocat(2,itypi,itypj) + al3 = alphisocat(3,itypi,itypj) + al4 = alphisocat(4,itypi,itypj) + csig = (1.0d0 & + / dsqrt(sigiso1cat(itypi, itypj)**2.0d0 & + + sigiso2cat(itypi,itypj)**2.0d0)) +!c! + pis = sig0headcat(itypi,itypj) + eps_head = epsheadcat(itypi,itypj) + Rhead_sq = Rhead * Rhead +!c! R1 - distance between head of ith side chain and tail of jth sidechain +!c! R2 - distance between head of jth side chain and tail of ith sidechain + R1 = 0.0d0 + R2 = 0.0d0 + DO k = 1, 3 +!c! Calculate head-to-tail distances needed by Epol + R1=R1+(ctail(k,2)-chead(k,1))**2 + R2=R2+(chead(k,2)-ctail(k,1))**2 + END DO +!c! Pitagoras + R1 = dsqrt(R1) + 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! Coulomb electrostatic interaction + Ecl = (332.0d0 * Qij) / Rhead +!c! derivative of Ecl is Gcl... + dGCLdR = (-332.0d0 * Qij ) / Rhead_sq + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq)) + Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0) + debkap=debaykapcat(itypi,itypj) + Egb = -(332.0d0 * Qij *& + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb +! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out +!c! Derivative of Egb is Ggb... + dGGBdFGB = -(-332.0d0 * Qij * & + (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)& + -(332.0d0 * Qij *& + (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb + dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) + dGGBdR = dGGBdFGB * dFGBdR +!c!------------------------------------------------------------------- +!c! Fisocav - isotropic cavity creation term +!c! or "how much energy it costs to put charged head in water" + pom = Rhead * csig + top = al1 * (dsqrt(pom) + al2 * pom - al3) + bot = (1.0d0 + al4 * pom**12.0d0) + botsq = bot * bot + FisoCav = top / bot +! write (*,*) "Rhead = ",Rhead +! write (*,*) "csig = ",csig +! write (*,*) "pom = ",pom +! write (*,*) "al1 = ",al1 +! write (*,*) "al2 = ",al2 +! write (*,*) "al3 = ",al3 +! write (*,*) "al4 = ",al4 +! write (*,*) "top = ",top +! write (*,*) "bot = ",bot +!c! Derivative of Fisocav is GCV... + dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2) + dbot = 12.0d0 * al4 * pom ** 11.0d0 + dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig +!c!------------------------------------------------------------------- +!c! Epol +!c! Polarization energy - charged heads polarize hydrophobic "neck" + MomoFac1 = (1.0d0 - chi1 * sqom2) + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR1 = ( R1 * R1 ) / MomoFac1 + RR2 = ( R2 * R2 ) / MomoFac2 + ee1 = exp(-( RR1 / (4.0d0 * a12sq) )) + ee2 = exp(-( RR2 / (4.0d0 * a12sq) )) + fgb1 = sqrt( RR1 + a12sq * ee1 ) + fgb2 = sqrt( RR2 + a12sq * ee2 ) + epol = 332.0d0 * eps_inout_fac * ( & + (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 )) +!c! epol = 0.0d0 + dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)& + / (fgb1 ** 5.0d0) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)& + / (fgb2 ** 5.0d0) + dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )& + / ( 2.0d0 * fgb1 ) + dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )& + / ( 2.0d0 * fgb2 ) + dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))& + * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 ) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))& + * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 ) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +!c! dPOLdR1 = 0.0d0 + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 +!c! dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Elj +!c! Lennard-Jones 6-12 interaction between heads + 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! These things do the dRdX derivatives, that is +!c! allow us to change what we see from function that changes with +!c! distance to function that changes with LOCATION (of the interaction +!c! site) + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1) + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,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+nres) ) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j) + facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) + facd4 = dtailcat(2,itypi,itypj) * vbld_inv(j) + +!c! Now we add appropriate partial derivatives (one in each dimension) + DO k = 1, 3 + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + 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)) + gradpepcatx(k,i) = gradpepcatx(k,i) & + - dGCLdR * pom& + - dGGBdR * pom& + - dGCVdR * pom& + - dPOLdR1 * hawk& + - 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& +! + dGGBdR * pom+ dGCVdR * pom& +! + dPOLdR1 * (erhead_tail(k,1)& +! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))& +! + dPOLdR2 * condor + dGLJdR * pom + + gradpepcat(k,i) = gradpepcat(k,i) & + - dGCLdR * erhead(k)& + - dGGBdR * erhead(k)& + - dGCVdR * erhead(k)& + - dPOLdR1 * erhead_tail(k,1)& + - dPOLdR2 * erhead_tail(k,2)& + - dGLJdR * erhead(k) + + gradpepcat(k,j) = gradpepcat(k,j) & + + dGCLdR * erhead(k) & + + dGGBdR * erhead(k) & + + dGCVdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dPOLdR2 * erhead_tail(k,2)& + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE eqq_cat !c!------------------------------------------------------------------- SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) use comm_momo @@ -26079,26 +27176,93 @@ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) DO k = 1, 3 - hawk = (erhead_tail(k,1) + & - facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + + gvdwx(k,i) = gvdwx(k,i) & + - dPOLdR1 * hawk + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR1 * (erhead_tail(k,1) & + -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) + + gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1) + gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1) + + END DO + RETURN + END SUBROUTINE eqn + SUBROUTINE enq(Epol) + use calc_data + use comm_momo + double precision facd3, adler,epol + alphapol2 = alphapol(itypj,itypi) +!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 Polarization energy + 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! Return the results +!c! (See comments in Eqq) + DO k = 1, 3 + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(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+nres))) gvdwx(k,i) = gvdwx(k,i) & - - dPOLdR1 * hawk - gvdwx(k,j) = gvdwx(k,j) & - + dPOLdR1 * (erhead_tail(k,1) & - -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) + - dPOLdR2 * (erhead_tail(k,2) & + -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR2 * condor - gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1) - gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1) + gvdwc(k,i) = gvdwc(k,i) & + - dPOLdR2 * erhead_tail(k,2) + gvdwc(k,j) = gvdwc(k,j) & + + dPOLdR2 * erhead_tail(k,2) END DO - RETURN - END SUBROUTINE eqn - SUBROUTINE enq(Epol) + RETURN + END SUBROUTINE enq + + SUBROUTINE enq_cat(Epol) use calc_data use comm_momo double precision facd3, adler,epol - alphapol2 = alphapol(itypj,itypi) + alphapol2 = alphapolcat(itypj,itypi) !c! R2 - distance between head of jth side chain and tail of ith sidechain R2 = 0.0d0 DO k = 1, 3 @@ -26132,34 +27296,36 @@ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 dPOLdOM2 = 0.0d0 + !c!------------------------------------------------------------------- !c! Return the results !c! (See comments in Eqq) DO k = 1, 3 erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) END DO - eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j) ) adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) facd2 = d2 * vbld_inv(j+nres) - 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+nres))) + + 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 RETURN - END SUBROUTINE enq + END SUBROUTINE enq_cat + SUBROUTINE eqd(Ecl,Elj,Epol) use calc_data use comm_momo @@ -26251,38 +27417,273 @@ - dPOLdR1 * hawk & - dGLJdR * pom - pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) - gvdwx(k,j) = gvdwx(k,j) & - + dGCLdR * pom & - + dPOLdR1 * (erhead_tail(k,1) & - -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) & - + dGLJdR * pom + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx(k,j) = gvdwx(k,j) & + + dGCLdR * pom & + + dPOLdR1 * (erhead_tail(k,1) & + -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) & + + dGLJdR * pom + + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR1 * erhead_tail(k,1) & + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE eqd + SUBROUTINE edq(Ecl,Elj,Epol) +! IMPLICIT NONE + use comm_momo + use calc_data + + double precision facd3, adler,ecl,elj,epol + alphapol2 = alphapol(itypj,itypi) + w1 = wqdip(1,itypi,itypj) + w2 = wqdip(2,itypi,itypj) + pis = sig0head(itypi,itypj) + eps_head = epshead(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) + 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 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+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) + eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(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+nres))) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx(k,i) = gvdwx(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+nres)) + gvdwx(k,j) = gvdwx(k,j) & + + dGCLdR * pom & + + dPOLdR2 * condor & + + dGLJdR * pom + + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR2 * erhead_tail(k,2) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE edq + + SUBROUTINE edq_cat(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 + 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!------------------------------------------------------------------- +!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+nres) ) + 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+nres) ) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j) + 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)) + 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 - gvdwc(k,i) = gvdwc(k,i) & - - dGCLdR * erhead(k) & - - dPOLdR1 * erhead_tail(k,1) & - - dGLJdR * erhead(k) + gradpepcat(k,i) = gradpepcat(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k) - gvdwc(k,j) = gvdwc(k,j) & - + dGCLdR * erhead(k) & - + dPOLdR1 * erhead_tail(k,1) & - + 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 eqd - SUBROUTINE edq(Ecl,Elj,Epol) -! IMPLICIT NONE - use comm_momo + 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 = alphapol(itypj,itypi) - w1 = wqdip(1,itypi,itypj) - w2 = wqdip(2,itypi,itypj) - pis = sig0head(itypi,itypj) - eps_head = epshead(itypi,itypj) + 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 @@ -26301,8 +27702,10 @@ !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) +! print *,"CO2", itypi,itypj +! print *,"CO?!.", w1,w2,Qj,om1 ECL = sparrow / Rhead**2.0d0 & - hawk / Rhead**4.0d0 !c!------------------------------------------------------------------- @@ -26311,9 +27714,10 @@ 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 !c Epol @@ -26344,50 +27748,57 @@ * (((-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+nres) ) - erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) - eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) - adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) - facd1 = d1 * vbld_inv(i+nres) - facd2 = d2 * vbld_inv(j+nres) - facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + 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+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) & - - dGCLdR * pom & - - dPOLdR2 * (erhead_tail(k,2) & - -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) & - - dGLJdR * pom + 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+nres)) - gvdwx(k,j) = gvdwx(k,j) & - + dGCLdR * pom & - + dPOLdR2 * condor & - + 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 - gvdwc(k,i) = gvdwc(k,i) & + gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( & - dGCLdR * erhead(k) & - dPOLdR2 * erhead_tail(k,2) & - - dGLJdR * erhead(k) + - dGLJdR * erhead(k)) + gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( & + - 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) END DO RETURN - END SUBROUTINE edq + END SUBROUTINE edq_cat_pep + SUBROUTINE edd(ECL) ! IMPLICIT NONE use comm_momo @@ -26493,51 +27904,270 @@ 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 ) +!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 ) +!c! distance from center of chain(?) to polar/charged head +!c! write (*,*) "istate = ", 1 +!c! write (*,*) "ii = ", 1 +!c! write (*,*) "jj = ", 1 + d1 = dhead(1, 1, itypi, itypj) + d2 = dhead(2, 1, itypi, itypj) +!c! ai*aj from Fgb + a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) +!c! a12sq = a12sq * a12sq +!c! charge of amino acid itypi is... + Qi = icharge(itypi) + Qj = icharge(itypj) + Qij = Qi * Qj +!c! chis1,2,12 + chis1 = chis(itypi,itypj) + chis2 = chis(itypj,itypi) + chis12 = chis1 * chis2 + sig1 = sigmap1(itypi,itypj) + sig2 = sigmap2(itypi,itypj) +!c! write (*,*) "sig1 = ", sig1 +!c! write (*,*) "sig2 = ", sig2 +!c! alpha factors from Fcav/Gcav + b1cav = alphasur(1,itypi,itypj) +! b1cav=0.0 + b2cav = alphasur(2,itypi,itypj) + b3cav = alphasur(3,itypi,itypj) + b4cav = alphasur(4,itypi,itypj) + wqd = wquad(itypi, itypj) +!c! used by Fgb + 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!------------------------------------------------------------------- +!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)!-dtailcat(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 = 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) +!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_cat + + 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 = 10 + 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 = 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 -!c! write (*,*) "istate = ", 1 -!c! write (*,*) "ii = ", 1 -!c! write (*,*) "jj = ", 1 - d1 = dhead(1, 1, itypi, itypj) - d2 = dhead(2, 1, itypi, itypj) + d1 = dheadcat(1, 1, itypi, itypj) + d2 = dheadcat(2, 1, itypi, itypj) !c! ai*aj from Fgb - a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) + a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj) !c! a12sq = a12sq * a12sq !c! charge of amino acid itypi is... - Qi = icharge(itypi) - Qj = icharge(itypj) - Qij = Qi * Qj + Qi = 0 + Qj = ichargecat(itypj) +! Qij = Qi * Qj !c! chis1,2,12 - chis1 = chis(itypi,itypj) - chis2 = chis(itypj,itypi) - chis12 = chis1 * chis2 - sig1 = sigmap1(itypi,itypj) - sig2 = sigmap2(itypi,itypj) -!c! write (*,*) "sig1 = ", sig1 -!c! write (*,*) "sig2 = ", sig2 + 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 = alphasur(1,itypi,itypj) -! b1cav=0.0 - b2cav = alphasur(2,itypi,itypj) - b3cav = alphasur(3,itypi,itypj) - b4cav = alphasur(4,itypi,itypj) - wqd = wquad(itypi, itypj) + 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 = epsintab(itypi,itypj) + 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)+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) @@ -26552,15 +28182,15 @@ !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,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) @@ -26585,7 +28215,7 @@ dPOLdOM1 = 0.0d0 dPOLdOM2 = 0.0d0 RETURN - END SUBROUTINE elgrad_init + END SUBROUTINE elgrad_init_cat_pep double precision function tschebyshev(m,n,x,y) implicit none @@ -26628,7 +28258,444 @@ 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 +!----------------------------------------------------------------------------- +!-----------------------------------------------------------------------------