X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=52fb0ff35b4211deb75f46c614fb14440f81b3cf;hb=93156c3542905e5b694a64496717b10b9c6a49a8;hp=4e043fe8faf529d49233e474f9bc01ab41e055a3;hpb=817262c16c50ae5848d85d1162e4977c700e1b6f;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 4e043fe..52fb0ff 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 @@ -383,8 +383,8 @@ time_Bcastw=time_Bcastw+MPI_Wtime()-time00 ! call chainbuild_cart endif -! print *,'Processor',myrank,' calling etotal ipot=',ipot -! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct + print *,'Processor',myrank,' calling etotal ipot=',ipot + print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else ! if (modecalc.eq.12.or.modecalc.eq.14) then ! call int_from_cart1(.false.) @@ -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 @@ -2842,6 +2848,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) @@ -3025,17 +3032,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(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-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 +4345,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 +4363,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 +4378,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 +4392,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 @@ -5040,9 +5054,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)) @@ -5073,7 +5087,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 @@ -5173,7 +5187,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) & @@ -11194,6 +11208,7 @@ #ifdef TIMING time01=MPI_Wtime() #endif +!#define DEBUG #ifdef DEBUG write (iout,*) "sum_gradient gvdwc, gvdwx" do i=1,nres @@ -11227,18 +11242,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 +12764,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 +22669,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 +22684,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 +22703,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 +22793,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 +22833,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 +22878,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 +22926,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 +23103,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)