! print *,"Processor",myrank," BROADCAST iorder"
! FG master sets up the WEIGHTS_ array which will be broadcast to the
! FG slaves as WEIGHTS array.
- ! weights_(1)=wsc
+ weights_(1)=wsc
weights_(2)=wscp
weights_(3)=welec
weights_(4)=wcorr
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
-! write(iout,*),"just befor eelec call"
+! print *,"just befor eelec call"
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-! write (iout,*) "ELEC calc"
+! print *, "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
eespp=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+! print *,"before ecatcat"
if (nfgtasks.gt.1) then
if (fg_rank.eq.0) then
call ecatcat(ecationcation)
! to calculate the el-loc multibody terms of various order.
!
!AL el mu=0.0d0
+
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti = nloctyp
+ else
iti = itype2loc(itype(i-2,1))
+ endif
else
iti=nloctyp
endif
!d write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
then
- call matmat2(CC(1,1,i-2),Ugder(1,1,i-2),CUgder(1,1,i-2))
+ call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
- call matvec2(Ctilde(1,1,i-2),obrot_der(1,i-2),Ctobrder(1,i-2))
+ 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(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-2),Ub2(1,i-2),CUgb2(1,i-2))
- call matvec2(CC(1,1,i-2),Ub2der(1,i-2),CUgb2der(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
+a23*gmuij1(2)&
+a32*gmuij1(3)&
+a33*gmuij1(4))&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
!c write(iout,*) "derivative over thatai"
!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
!c & a33*gmuij1(4)
+a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
geel_loc_ij*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
!c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)&
gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
geel_loc_ji*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
geel_loc_ji=&
+a22*gmuji2(1)&
!c & a33*gmuji2(4)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
geel_loc_ji*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
#endif
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
#ifdef TIMING
time01=MPI_Wtime()
#endif
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
do i=1,nres
! write (iout,'(i5,3f10.5,2x,f10.5)')
! & i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
! enddo
- write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
- (gvdwc_scpp(j,i),j=1,3)
- enddo
- write (iout,*) "gelc_long gvdwpp gel_loc_long"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
- (gelc_loc_long(j,i),j=1,3)
- enddo
+! write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+! (gvdwc_scpp(j,i),j=1,3)
+! enddo
+! write (iout,*) "gelc_long gvdwpp gel_loc_long"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+! (gelc_loc_long(j,i),j=1,3)
+! enddo
call flush(iout)
#endif
#ifdef SPLITELE
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=2.0D-5
+ aincr=1.0D-7
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
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
+ 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, &
ecation_prot=0.0d0
! first lets calculate interaction with peptide groups
if (nres_molec(5).eq.0) return
- wconst=78
- 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
do i=1,4
itmp=itmp+nres_molec(i)
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)
E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
2*wvan2**6*Irsistp)
ecation_prot = ecation_prot+E1+E2
+! print *,"ecatprot",i,j,ecation_prot,rcpm
dE1dr = -2*costhet*wdip*Irthrp-&
(4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
dE2dr = 3*wquad1*wquad2*Irfourp- &
enddo
cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
do j=itmp+1,itmp+nres_molec(5)
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
endif
! enddo
! enddo
- if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+ if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+ ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+ (itype(i,1).eq.25))) then
if(itype(i,1).eq.16) then
inum=1
else
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of calcium with ASP/GLU
- wh2o=78
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ ndivi=0.5
+ else
+ ndivi=1.0
+ endif
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndivi*ndiv
wc = vcatprm(1)
wc=wc/wh2o
wdip =vcatprm(2)
v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of ASN/GLN with calcium
- wh2o=78
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndiv
wdip =vcatprm(2)
wdip=wdip/wh2o
wquad1 =vcatprm(3)