From: Adam Sieradzan Date: Wed, 16 Aug 2017 06:50:06 +0000 (+0200) Subject: working gradient X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=340832b9074a6903554ccfcf3cac2abf65e09d3a;p=unres4.git working gradient --- diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index 34f8a2a..8e6c0af 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -300,6 +300,7 @@ ncont=0 ees=0.0 evdw=0.0 + print *, "nntt,nct",nnt,nct-2 do 1 i=nnt,nct-2 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1 xi=c(1,i) @@ -492,6 +493,8 @@ enddo call elecont(lprint,ncont,icont) + print *,"after elecont" + if (nres_molec(1).eq.0) return ! finding parallel beta !d write (iout,*) '------- looking for parallel beta -----------' @@ -742,7 +745,7 @@ write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'" write(12,'(a20)') "XMacStand ribbon.mac" - + if (nres_molec(1).eq.0) return write(iout,*) 'UNRES seq:' do j=1,nbfrag write(iout,*) 'beta ',(bfrag(i,j),i=1,4) @@ -815,6 +818,7 @@ ! & obr,non_conv) ! rms=dsqrt(rms) call rmsd(rms) + print *,"before contact" !elte(iout,*) "rms_nacc before contact" call contact(.false.,ncont,icont,co) frac=contact_fract(ncont,ncont_ref,icont,icont_ref) @@ -877,6 +881,7 @@ ! else ! do kkk=1,nperm iatom=0 + print *,nz_start,nz_end,nstart_seq-nstart_sup do i=nz_start,nz_end iatom=iatom+1 iti=itype(i,1) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index d70d144..24ea5e5 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -728,7 +728,9 @@ ispp=4 !?? wham ispp=2 #ifdef MPI ! Now partition the electrostatic-interaction array - if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then + if (nres_molec(1).eq.0) then + npept=0 + elseif (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then npept=nres_molec(1)-nnt-1 else npept=nres_molec(1)-nnt @@ -787,7 +789,7 @@ ijunk,ielstart_nucl(i),ielend_nucl(i),*113) enddo ! i 113 continue - if (iatel_s.eq.0) iatel_s=1 + if (iatel_s_nucl.eq.0) iatel_s_nucl=1 nele_int_tot_vdw=(npept-2)*(npept-2+1)/2 ! write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw @@ -810,7 +812,7 @@ if (iatel_s_vdw.eq.0) iatel_s_vdw=1 15 continue if (iatel_s.eq.0) iatel_s=1 - + if (iatel_s_vdw.eq.0) iatel_s_vdw=1 nele_int_tot_vdw_nucl=(npept_nucl-2)*(npept_nucl-2+1)/2 ! write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw call int_bounds(nele_int_tot_vdw_nucl,my_ele_inds_vdw_nucl,& @@ -831,7 +833,7 @@ ! write (iout,*) i," ielstart_vdw",ielstart_vdw(i), ! & " ielend_vdw",ielend_vdw(i) enddo ! i - if (iatel_s_vdw.eq.0) iatel_s_vdw=1 + if (iatel_s_vdw.eq.0) iatel_s_vdw_nucl=1 115 continue #else @@ -934,6 +936,7 @@ enddo ! i 114 continue print *, "after inloop3",iatscp_s_nucl,iatscp_e_nucl + if (iatscp_s_nucl.eq.0) iatscp_s_nucl=1 #else iatscp_s=nnt iatscp_e=nct_molec(1)-1 diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index bbc2ed2..baf7f55 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -66,7 +66,7 @@ wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,& wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, & wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,& - wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpsb + wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb #ifdef CLUSTER real(kind=8) :: scalscp #endif diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 881497c..9d92c63 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -128,7 +128,8 @@ !-----------------------------NUCLEIC GRADIENT real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl, & gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,& - gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl + gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl,& + gvdwpp_nucl ! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2) real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,& gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres) @@ -274,6 +275,20 @@ weights_(17)=wbond weights_(18)=scal14 weights_(21)=wsccor + weights_(26)=wvdwpp_nucl + weights_(27)=welpp + weights_(28)=wvdwpsb + weights_(29)=welpsb + weights_(30)=wvdwsb + weights_(31)=welsb + weights_(32)=wbond_nucl + weights_(33)=wang_nucl + weights_(34)=wsbloc + weights_(35)=wtor_nucl + weights_(36)=wtor_d_nucl + weights_(37)=wcorr_nucl + weights_(38)=wcorr3_nucl + ! FG Master broadcasts the WEIGHTS_ array call MPI_Bcast(weights_(1),n_ene,& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) @@ -300,6 +315,20 @@ wbond=weights(17) scal14=weights(18) wsccor=weights(21) + wvdwpp_nucl =weights(26) + welpp =weights(27) + wvdwpsb=weights(28) + welpsb =weights(29) + wvdwsb =weights(30) + welsb =weights(31) + wbond_nucl =weights(32) + wang_nucl =weights(33) + wsbloc =weights(34) + wtor_nucl =weights(35) + wtor_d_nucl =weights(36) + wcorr_nucl =weights(37) + wcorr3_nucl =weights(38) + endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -365,7 +394,7 @@ if (shield_mode.eq.2) then call set_shield_fac2 endif - print *,"AFTER EGB",ipot,evdw +! print *,"AFTER EGB",ipot,evdw !mc !mc Sep-06: egb takes care of dynamic ss bonds too !mc @@ -431,7 +460,7 @@ ! Calculate the bond-stretching energy ! call ebond(estr) - print *,"EBOND",estr +! print *,"EBOND",estr ! write(iout,*) "in etotal afer ebond",ipot ! @@ -448,6 +477,7 @@ call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif ! print *,"Processor",myrank," computed UB" ! @@ -549,7 +579,7 @@ etube=0.0d0 endif !-------------------------------------------------------- - print *,"before",ees,evdw1,ecorr +! print *,"before",ees,evdw1,ecorr call ebond_nucl(estr_nucl) call ebend_nucl(ebe_nucl) call etor_nucl(etors_nucl) @@ -559,7 +589,7 @@ call esb(esbloc) call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1) - print *,"after ebend", ebe_nucl +! print *,"after ebend", ebe_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -753,7 +783,7 @@ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& +Eafmforce+ethetacnstr & +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& - +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& + +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl #else @@ -765,7 +795,7 @@ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& +Eafmforce+ethetacnstr & +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& - +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& + +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl #endif @@ -952,7 +982,7 @@ edihcnstr,ethetacnstr,ebr*nss,& Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, & - evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,& + 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, & @@ -1008,7 +1038,7 @@ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & etube,wtube, & estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,& - evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb& + 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, & @@ -3067,7 +3097,7 @@ ! ! 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 +! print *,"iatel_s,iatel_e,",iatel_s,iatel_e do i=iatel_s,iatel_e if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle dxi=dc(1,i) @@ -5461,6 +5491,8 @@ if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) enddo +! print *,ithetaconstr_start,ithetaconstr_end,"TU" + ! Ufff.... We've done all this!!! return end subroutine ebend @@ -5762,7 +5794,7 @@ !-----------thete constrains ! if (tor_mode.ne.2) then ethetacnstr=0.0d0 -!C print *,ithetaconstr_start,ithetaconstr_end,"TU" +! print *,ithetaconstr_start,ithetaconstr_end,"TU" do i=ithetaconstr_start,ithetaconstr_end itheta=itheta_constr(i) thetiii=theta(itheta) @@ -10410,7 +10442,12 @@ +wturn4*gshieldc_t4(j,i)& +wel_loc*gshieldc_ll(j,i)& +wtube*gg_tube(j,i) & - +wbond_nucl*gradb_nucl(j,i) + +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ & + wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ & + wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ & + wcorr_nucl*gradcorr_nucl(j,i)& + +wcorr3_nucl*gradcorr3_nucl(j,i) + enddo enddo #else @@ -10433,8 +10470,11 @@ +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i)& +wtube*gg_tube(j,i) & - +wbond_nucl*gradb_nucl(j,i) - + +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ & + wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ & + wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ & + wcorr_nucl*gradcorr_nucl(j,i) + +wcorr3_nucl*gradcorr3_nucl(j,i) enddo enddo #endif @@ -10591,7 +10631,46 @@ +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & +wtube*gg_tube(j,i) & - +wbond_nucl*gradb_nucl(j,i) + +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)& + +wvdwpsb*gvdwpsb1(j,i))& + +wbond_nucl*gradb_nucl(j,i)+wsbloc*gsbloc(j,i) + +! if ((i.le.2).and.(i.ge.1)) print *,gradc(j,i,icg),& +! gradbufc(j,i),welec*gelc(j,i), & +! wel_loc*gel_loc(j,i), & +! wscp*gvdwc_scpp(j,i), & +! welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), & +! wel_loc*gel_loc_long(j,i), & +! wcorr*gradcorr_long(j,i), & +! wcorr5*gradcorr5_long(j,i), & +! wcorr6*gradcorr6_long(j,i), & +! wturn6*gcorr6_turn_long(j,i), & +! wbond*gradb(j,i), & +! wcorr*gradcorr(j,i), & +! wturn3*gcorr3_turn(j,i), & +! wturn4*gcorr4_turn(j,i), & +! wcorr5*gradcorr5(j,i), & +! wcorr6*gradcorr6(j,i), & +! wturn6*gcorr6_turn(j,i), & +! wsccor*gsccorc(j,i) & +! ,wscloc*gscloc(j,i) & +! ,wliptran*gliptranc(j,i) & +! ,gradafm(j,i) & +! ,welec*gshieldc(j,i) & +! ,welec*gshieldc_loc(j,i) & +! ,wcorr*gshieldc_ec(j,i) & +! ,wcorr*gshieldc_loc_ec(j,i) & +! ,wturn3*gshieldc_t3(j,i) & +! ,wturn3*gshieldc_loc_t3(j,i) & +! ,wturn4*gshieldc_t4(j,i) & +! ,wturn4*gshieldc_loc_t4(j,i) & +! ,wel_loc*gshieldc_ll(j,i) & +! ,wel_loc*gshieldc_loc_ll(j,i) & +! ,wtube*gg_tube(j,i) & +! ,wbond_nucl*gradb_nucl(j,i) & +! ,wvdwpp_nucl*gvdwpp_nucl(j,i),welpp*gelpp(j,i),& +! wvdwpsb*gvdwpsb1(j,i)& +! ,wbond_nucl*gradb_nucl(j,i),wsbloc*gsbloc(j,i) @@ -10627,7 +10706,10 @@ +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & +wtube*gg_tube(j,i) & - +wbond_nucl*gradb_nucl(j,i) + +wbond_nucl*gradb_nucl(j,i) & + +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)& + +wvdwpsb*gvdwpsb1(j,i))& + +wsbloc*gsbloc(j,i) @@ -10645,10 +10727,12 @@ +wturn4*gshieldx_t4(j,i) & +wel_loc*gshieldx_ll(j,i)& +wtube*gg_tube_sc(j,i) & - +wbond_nucl*gradbx_nucl(j,i) - - - + +wbond_nucl*gradbx_nucl(j,i) & + +wvdwsb*gvdwsbx(j,i) & + +welsb*gelsbx(j,i) & + +wcorr_nucl*gradxorr_nucl(j,i)& + +wcorr3_nucl*gradxorr3_nucl(j,i) & + +wsbloc*gsblocx(j,i) enddo enddo #ifdef DEBUG @@ -10675,7 +10759,7 @@ #ifdef MPI if (nfgtasks.gt.1) then do j=1,3 - do i=1,nres + do i=0,nres gradbufc(j,i)=gradc(j,i,icg) gradbufx(j,i)=gradx(j,i,icg) enddo @@ -10704,7 +10788,7 @@ time00=MPI_Wtime() call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) - call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,& + call MPI_Reduce(gradbufx(1,0),gradx(1,0,icg),3*nres+3,& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) @@ -10713,6 +10797,7 @@ 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 @@ -11650,7 +11735,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo call zerograd call etotal_short(energia) -!el call enerprint(energia) + call enerprint(energia) call flush(iout) write (iout,*) "enter cartgrad" call flush(iout) @@ -11857,6 +11942,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) +! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo @@ -11934,6 +12020,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (.not.split_ene) then call etotal(energia1) etot1=energia1(0) +! call enerprint(energia1) else !- split gradient call etotal_long(energia1) @@ -12062,11 +12149,11 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call var_to_geom(nvar,x) call chainbuild icall=1 - print *,'ICG=',ICG +! print *,'ICG=',ICG call etotal(energia) etot = energia(0) !el call enerprint(energia) - print *,'ICG=',ICG +! print *,'ICG=',ICG #ifdef MPL if (MyID.ne.BossID) then call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID) @@ -16070,6 +16157,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do j=1,3 gcart(j,i)=gradc(j,i,icg) gxcart(j,i)=gradx(j,i,icg) +! if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg) enddo #ifdef DEBUG write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),& @@ -16080,6 +16168,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' time01=MPI_Wtime() #endif call int_to_cart +! print *,"gcart_two",gcart(2,2),gradc(2,2,icg) + #ifdef TIMING time_inttocart=time_inttocart+MPI_Wtime()-time01 #endif @@ -16230,6 +16320,25 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gradafm(j,i)=0.0d0 gradb_nucl(j,i)=0.0d0 gradbx_nucl(j,i)=0.0d0 + gvdwpp_nucl(j,i)=0.0d0 + gvdwpp(j,i)=0.0d0 + gelpp(j,i)=0.0d0 + gvdwpsb(j,i)=0.0d0 + gvdwpsb1(j,i)=0.0d0 + gvdwsbc(j,i)=0.0d0 + gvdwsbx(j,i)=0.0d0 + gelsbc(j,i)=0.0d0 + gradcorr_nucl(j,i)=0.0d0 + gradcorr3_nucl(j,i)=0.0d0 + gradxorr_nucl(j,i)=0.0d0 + gradxorr3_nucl(j,i)=0.0d0 + gelsbx(j,i)=0.0d0 + gsbloc(j,i)=0.0d0 + gsblocx(j,i)=0.0d0 + enddo + enddo + do i=0,nres + do j=1,3 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo @@ -19604,6 +19713,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(gradxorr_nucl(3,-1:nres)) allocate(gradcorr3_nucl(3,-1:nres)) allocate(gradxorr3_nucl(3,-1:nres)) + allocate(gvdwpp_nucl(3,-1:nres)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) @@ -19813,14 +19923,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),& vbldp0_nucl,diff,AKP_nucl*diff*diff estr_nucl=estr_nucl+diff*diff - print *,estr_nucl +! print *,estr_nucl do j=1,3 gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i) enddo !c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) enddo estr_nucl=0.5d0*AKP_nucl*estr_nucl - print *,"partial sum", estr_nucl,AKP_nucl +! print *,"partial sum", estr_nucl,AKP_nucl if (energy_dec) & write (iout,*) "ibondp_start,ibondp_end",& @@ -19839,7 +19949,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, & AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff - print *,estr_nucl +! print *,estr_nucl do j=1,3 gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) enddo @@ -19869,7 +19979,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo estr_nucl=estr_nucl+uprod/usum do j=1,3 - gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) + gradbx_nucl(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) enddo endif enddo @@ -19882,7 +19992,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble - logical :: lprn=.true., lprn1=.false. + logical :: lprn=.false., lprn1=.false. !el local variables integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai @@ -20166,7 +20276,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !c !c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 !c - print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl +! print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl do i=iatel_s_nucl,iatel_e_nucl if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle dxi=dc(1,i) @@ -20253,8 +20363,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ggg(2)=fac*yj ggg(3)=fac*zj do k=1,3 - gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) - gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) + gvdwpp_nucl(k,i)=gvdwpp_nucl(k,i)-ggg(k) + gvdwpp_nucl(k,j)=gvdwpp_nucl(k,j)+ggg(k) enddo !c phoshate-phosphate electrostatic interactions rij=dsqrt(rij) @@ -20282,7 +20392,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=nnt,nct !c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3) do k=1,3 - gvdwpp(k,i)=6*gvdwpp(k,i) + gvdwpp_nucl(k,i)=6*gvdwpp_nucl(k,i) !c gelpp(k,i)=332.0d0*gelpp(k,i) gelpp(k,i)=AEES*gelpp(k,i) enddo @@ -20311,7 +20421,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e eelpsb=0.0d0 evdwpsb=0.0d0 - print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl +! print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl do i=iatscp_s_nucl,iatscp_e_nucl if (itype(i,2).eq.ntyp1_molec(2) & .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle @@ -20910,8 +21020,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo do j = 1,3 z_prime(j) = -uz(j,i-1) +! z_prime(j)=0.0 enddo - + xx=0.0d0 yy=0.0d0 zz=0.0d0 @@ -20941,8 +21052,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #endif sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) esbloc = esbloc + sumene - if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz - if (energy_dec) write(iout,*) "x",(x(k),k=1,9) + sumene2= enesc_nucl(x,xx,yy,0.0d0,cost2tab(i+1),sint2tab(i+1)) +! print *,"enecomp",sumene,sumene2 +! if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz +! if (energy_dec) write(iout,*) "x",(x(k),k=1,9) #ifdef DEBUG write (2,*) "x",(x(k),k=1,9) !C @@ -21054,12 +21167,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !c & dt_dci(k) !c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ", !c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k) - gsbloc(k,i-1)=gsbloc(k,i-1)+de_dxx*dxx_ci1(k) & - +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k) - gsbloc(k,i)=gsbloc(k,i)+de_dxx*dxx_Ci(k) & - +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k) + gsbloc(k,i-1)=gsbloc(k,i-1)+(de_dxx*dxx_ci1(k) & + +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)) + gsbloc(k,i)=gsbloc(k,i)+(de_dxx*dxx_Ci(k) & + +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)) gsblocx(k,i)= de_dxx*dxx_XYZ(k)& +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k) +! print *,i,de_dxx*dxx_ci1(k)+de_dyy*dyy_ci1(k),de_dzz*dzz_ci1(k)*2 enddo !c write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3), !c & (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3) @@ -21233,21 +21347,21 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ecorr=0.0D0 ecorr3=0.0d0 !C Remove the loop below after debugging !!! - do i=nnt_molec(2),nct_molec(2) - do j=1,3 - gradcorr_nucl(j,i)=0.0D0 - gradxorr_nucl(j,i)=0.0D0 - gradcorr3_nucl(j,i)=0.0D0 - gradxorr3_nucl(j,i)=0.0D0 - enddo - enddo - print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl +! do i=nnt_molec(2),nct_molec(2) +! do j=1,3 +! gradcorr_nucl(j,i)=0.0D0 +! gradxorr_nucl(j,i)=0.0D0 +! gradcorr3_nucl(j,i)=0.0D0 +! gradxorr3_nucl(j,i)=0.0D0 +! enddo +! enddo +! print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl !C Calculate the local-electrostatic correlation terms do i=iatsc_s_nucl,iatsc_e_nucl i1=i+1 num_conti=num_cont_hb(i) num_conti1=num_cont_hb(i+1) - print *,i,num_conti,num_conti1 +! print *,i,num_conti,num_conti1 do jj=1,num_conti j=jcont_hb(jj,i) do kk=1,num_conti1 diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index e0b37d3..653e595 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -3240,6 +3240,8 @@ ! calculating dE/ddc1 !el local variables integer :: j,i +! print *,"gloc",gloc(:,:) +! print *, "gcart",gcart(:,:) if (nres.lt.3) go to 18 do j=1,3 gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) & @@ -3331,7 +3333,8 @@ enddo ! The side-chain vector derivatives do i=2,nres-1 - if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if(itype(i,1).ne.10 .and. & + itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then do j=1,3 gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) & +gloc(ialph(i,1)+nside,icg)*domega(j,3,i) @@ -3350,7 +3353,8 @@ ! enddo if (nres.lt.2) return if ((nres.lt.3).and.(itype(1,1).eq.10)) return - if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then + if ((itype(1,1).ne.10).and. & + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then do j=1,3 !c Derviative was calculated for oposite vector of side chain therefore ! there is "-" sign before gloc_sc @@ -3358,7 +3362,8 @@ dtauangle(j,1,1,3) gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* & dtauangle(j,1,2,3) - if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then + if ((itype(2,1).ne.10).and. & + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then gxcart(j,1)= gxcart(j,1) & -gloc_sc(3,0,icg)*dtauangle(j,3,1,3) gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* & @@ -3366,7 +3371,8 @@ endif enddo endif - if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) & + if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.& + (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) & then do j=1,3 gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4) @@ -3378,10 +3384,11 @@ ! Calculating the remainder of dE/ddc2 do j=1,3 - if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then + if((itype(2,1).ne.10).and. & + (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ & gloc_sc(3,0,icg)*dtauangle(j,3,3,3) - if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) & + if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) & then gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4) !c the - above is due to different vector direction @@ -3395,7 +3402,8 @@ ! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx" endif endif - if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then + if ((itype(1,1).ne.10).and.& + (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3) ! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3) endif @@ -3415,7 +3423,8 @@ do i=3,nres-2 do j=1,3 ! write(iout,*) "before", gcart(j,i) - if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then + if ((itype(i,1).ne.10).and.& + (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) & *dtauangle(j,2,3,i+1) & -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2) @@ -3454,7 +3463,8 @@ ! Setting dE/ddnres-1 if(nres.ge.4) then do j=1,3 - if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then + if ((itype(nres-1,1).ne.10).and.& + (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) & *dtauangle(j,2,3,nres) ! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg), @@ -3463,18 +3473,20 @@ gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) & *dtauangle(j,3,3,nres) endif - if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then + if ((itype(nres,1).ne.10).and.& + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) & *dtauangle(j,3,1,nres+1) gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) & *dtauangle(j,3,2,nres+1) endif endif - if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then + if ((itype(nres-2,1).ne.10).and.& + (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* & dtauangle(j,1,3,nres) endif - if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then + if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* & dtauangle(j,2,2,nres+1) ! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg), @@ -3484,7 +3496,7 @@ endif ! Settind dE/ddnres if ((nres.ge.3).and.(itype(nres,1).ne.10).and. & - (itype(nres,1).ne.ntyp1))then + (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then do j=1,3 gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) & *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) & @@ -3492,6 +3504,7 @@ enddo endif ! The side-chain vector derivatives +! print *,"gcart",gcart(:,:) return end subroutine int_to_cart #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) diff --git a/source/unres/io.f90 b/source/unres/io.f90 index ee588d0..a86fc0b 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -720,7 +720,7 @@ ! integer :: rescode ! double precision x(maxvar) character(len=256) :: pdbfile - character(len=480) :: weightcard + character(len=800) :: weightcard character(len=80) :: weightcard_t!,ucase ! integer,dimension(:),allocatable :: itype_pdb !(maxres) ! common /pizda/ itype_pdb @@ -769,6 +769,19 @@ call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) + call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,1.0D0) + call reada(weightcard,'WELPP',welpp,1.0d0) + call reada(weightcard,'WVDWPSB',wvdwpsb,1.0d0) + call reada(weightcard,'WELPSB',welpsb,1.0D0) + call reada(weightcard,'WVDWSB',wvdwsb,1.0d0) + call reada(weightcard,'WELSB',welsb,1.0D0) + call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0) + call reada(weightcard,'WANG_NUCL',wang_nucl,1.0D0) + call reada(weightcard,'WSBLOC',wsbloc,1.0D0) + call reada(weightcard,'WTOR_NUCL',wtor_nucl,1.0D0) + call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,1.0D0) + call reada(weightcard,'WCORR_NUCL',wcorr_nucl,1.0D0) + call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) @@ -1188,8 +1201,8 @@ #endif nct=nres print *,'NNT=',NNT,' NCT=',NCT - if (itype(1,1).eq.ntyp1) nnt=2 - if (itype(nres,1).eq.ntyp1) nct=nct-1 + if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2 + if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1 if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup diff --git a/source/unres/unres.f90 b/source/unres/unres.f90 index 40c0140..e84e1d0 100644 --- a/source/unres/unres.f90 +++ b/source/unres/unres.f90 @@ -346,7 +346,7 @@ write (iout,'(a,i20)') '# of energy evaluations:',nfun+1 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals else - print *,'refstr=',refstr + print *,'refstr=',refstr,frac,frac_nn,co if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) print *,"after rms_nac_ncc" call briefout(0,etot)