X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=3870bb0aaeef6d24551f54ef0c7b9693d9426464;hb=da0ffdc0802840f3e36a2f41166ad72ced8f7845;hp=4e043fe8faf529d49233e474f9bc01ab41e055a3;hpb=8cd7afe8e87106d4cdd4b6361171845086971607;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 4e043fe..3870bb0 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -294,7 +294,7 @@ ! 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 @@ -610,9 +610,9 @@ .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 @@ -831,6 +831,7 @@ 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) @@ -2727,13 +2728,18 @@ ! 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 @@ -3025,17 +3031,17 @@ !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)) @@ -4338,7 +4344,9 @@ +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) @@ -4354,7 +4362,9 @@ +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)& @@ -4367,7 +4377,9 @@ 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)& @@ -4379,7 +4391,8 @@ !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 @@ -11194,6 +11207,7 @@ #ifdef TIMING time01=MPI_Wtime() #endif +!#define DEBUG #ifdef DEBUG write (iout,*) "sum_gradient gvdwc, gvdwx" do i=1,nres @@ -11227,18 +11241,18 @@ ! 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 @@ -12749,7 +12763,7 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=2.0D-5 + aincr=1.0D-7 write(iout,*) 'Calling CHECK_ECARTINT.',aincr nf=0 icall=0 @@ -22654,7 +22668,8 @@ 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, & @@ -22668,17 +22683,6 @@ 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) @@ -22698,6 +22702,23 @@ 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) @@ -22771,6 +22792,7 @@ 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- & @@ -22810,6 +22832,9 @@ 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) @@ -22852,7 +22877,10 @@ 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 @@ -22897,7 +22925,15 @@ ! 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) @@ -23066,7 +23102,10 @@ 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)