X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=badacd2fb7c07a79296eb61c383dd59001cc8963;hb=df8ac3297a3b0c878c3da92c9d3640f57b76cfd0;hp=205f78a063dafb7ffcffd877c463e97e5151ad67;hpb=d24fd92f8762d3b48a94c699ef90e0ab37d3c204;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 205f78a..badacd2 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -195,8 +195,8 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) - else + call ebend(ebe) + else ebe=0 endif c print *,"Processor",myrank," computed UB" @@ -215,7 +215,7 @@ cd print *,'nterm=',nterm else etors=0 edihcnstr=0 - endif + endif if (constr_homology.ge.1) then call e_modeller(ehomology_constr) @@ -233,9 +233,9 @@ C C 6/23/01 Calculate double-torsional energy C if (wtor_d.gt.0) then - call etor_d(etors_d) + call etor_d(etors_d) else - etors_d=0 + etors_d=0 endif c print *,"Processor",myrank," computed Utord" C @@ -277,14 +277,14 @@ c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr else Esaxs_constr = 0.0d0 - endif + endif c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy C after the equilibration time if(usampl.and.totT.gt.eq_time) then call EconstrQ - call Econstr_back + call Econstr_back else Uconst=0.0d0 Uconst_back=0.0d0 @@ -527,7 +527,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef DEBUG write (iout,*) "sum_gradient gsaxsc, gsaxsx" do i=0,nres - write (iout,'(i3,3e15.5,5x,3e15.5)') + write (iout,'(i3,3e15.5,5x,3e15.5)') & i,(gsaxsc(j,i),j=1,3),(gsaxsx(j,i),j=1,3) enddo call flush(iout) @@ -738,7 +738,7 @@ c enddo gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ - & welec*gelc_long(j,i) + + & welec*gelc_long(j,i)+ & wel_loc*gel_loc_long(j,i)+ & wcorr*gcorr_long(j,i)+ & wcorr5*gradcorr5_long(j,i)+ @@ -1530,6 +1530,7 @@ C include 'COMMON.SBRIDGE' logical lprn integer xshift,yshift,zshift + evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -2779,7 +2780,7 @@ c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) - enddo + enddo #else if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) @@ -2796,10 +2797,10 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) b2(2,i-2)=b(4,iti) - b1tilde(1,i-2)=b1(1,i-2) - b1tilde(2,i-2)=-b1(2,i-2) - b2tilde(1,i-2)=b2(1,i-2) - b2tilde(2,i-2)=-b2(2,i-2) + b1tilde(1,i-2)= b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)= b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) EE(1,2,i-2)=eeold(1,2,iti) EE(2,1,i-2)=eeold(2,1,iti) EE(2,2,i-2)=eeold(2,2,iti) @@ -2893,7 +2894,7 @@ cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) cd write (iout,*) 'b2',b2(:,iti) cd write (iout,*) "phi(",i,")=",phi(i)," sin1",sin1," cos1",cos1 -cd write (iout,*) 'Ug',Ug(:,:,i-2) +cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) @@ -3730,6 +3731,7 @@ C erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * @@ -5189,10 +5191,16 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr -cd write(iout,*)'link_start=',link_start,' link_end=',link_end +c write(iout,*)'link_start=',link_start,' link_end=',link_end, +c & " constr_dist",constr_dist if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@ -5217,38 +5225,75 @@ cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then C 15/02/13 CC dynamic SSbond - additional check - if (ii.gt.nres - & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then - call ssbond_ene(iii,jjj,eij) - ehpb=ehpb+2*eij + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then + call ssbond_ene(iii,jjj,eij) + ehpb=ehpb+2*eij endif cd write (iout,*) "eij",eij - else +cd & ' waga=',waga,' fac=',fac +! else if (ii.gt.nres .and. jj.gt.nres) then + else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (irestr_type(i).eq.11) then + ehpb=ehpb+fordepth(i)!**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)!**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)') + & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), + & ehpb,irestr_type(i) + else if (irestr_type(i).eq.10) then +c AL 6//19/2018 cross-link restraints + xdis = 0.5d0*(dd/forcon(i))**2 + expdis = dexp(-xdis) +c aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i) + aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i) +c write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux, +c & " wboltzd",wboltzd + ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux) +c fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i)) + fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i)) + & *expdis/(aux*forcon(i)**2) + if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)') + & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), + & -wboltzd*xlscore(i)*dlog(aux),irestr_type(i) + else if (irestr_type(i).eq.2) then +c Quartic restraints + ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') + & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), + & forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(i) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd + else +c Quadratic restraints rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+0.5d0*waga*rdis*rdis + if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') + & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), + & 0.5d0*waga*rdis*rdis,irestr_type(i) C C Evaluate gradient. C fac=waga*rdis/dd -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo + endif +c Calculate Cartesian gradient + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). if (iii.lt.ii) then - do j=1,3 - ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) - ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) - enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo endif cgrad do j=iii,jjj-1 cgrad do k=1,3 @@ -5261,7 +5306,6 @@ cgrad enddo enddo endif enddo - ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -5400,6 +5444,7 @@ C NO vbldp0 is the equlibrium lenght of spring for peptide group c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) c endif enddo + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included @@ -5580,6 +5625,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. 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)+gloc(nphi+i-2,icg) enddo + C Ufff.... We've done all this!!! return end @@ -5896,6 +5942,7 @@ c lprn1=.false. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) enddo + return end #endif @@ -6895,7 +6942,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6914,9 +6961,13 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d j = jres_homo(ii) dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij + nexl=0 do k=1,constr_homology c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii) - if(.not.l_homo(k,ii)) cycle + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -6955,7 +7006,15 @@ c write (iout,* )"min_odl",min_odl write (iout,*) "distancek",(distancek(k),k=1,constr_homology) write (iout,* )"min_odl",min_odl #endif +#ifdef OLDRESTR odleg2=0.0d0 +#else + if (waga_dist.ge.0.0d0) then + odleg2=nexl + else + odleg2=0.0d0 + endif +#endif do k=1,constr_homology c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ @@ -7110,8 +7169,11 @@ c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= c & -(6.28318-dih_diff(i,k)) c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= c & 6.28318+dih_diff(i,k) - +#ifdef OLD_DIHED kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#else + kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#endif c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) gdih(k)=dexp(kat3) kat2=kat2+gdih(k) @@ -7138,7 +7200,11 @@ c ---------------------------------------------------------------------- sum_gdih=kat2 sum_sgdih=0.0d0 do k=1,constr_homology +#ifdef OLD_DIHED sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd +#else + sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd +#endif c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle sum_sgdih=sum_sgdih+sgdih enddo @@ -7220,7 +7286,7 @@ c utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? gtheta(k)=dexp(utheta_i) ! + min_utheta_i? - gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk) + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) c Gradient for single Gaussian restraint in subr Econstr_back c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) c @@ -7292,7 +7358,9 @@ c c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? c uscdiffk(k)=usc_diff(i) guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff - guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk) +c write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i), +c & " guscdiff2",guscdiff2(k) + guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk) c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), c & xxref(j),yyref(j),zzref(j) enddo @@ -8566,7 +8634,7 @@ c & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees, c & 'gradcorr_long' C Calculate the multi-body contribution to energy. -c ecorr=ecorr+ekont*ees +C ecorr=ecorr+ekont*ees C Calculate multi-body contributions to the gradient. coeffpees0pij=coeffp*ees0pij coeffmees0mij=coeffm*ees0mij @@ -11050,6 +11118,42 @@ C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist end c---------------------------------------------------------------------------- + double precision function sscale2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) +c write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb +c write (2,*) "rr",rr + if(rr.lt.r_cut-rlamb) then + sscale2=1.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale2=0d0 + endif + return + end +C----------------------------------------------------------------------- + double precision function sscalgrad2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) + if(rr.lt.r_cut-rlamb) then + sscalgrad2=0.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + if (r.ge.r0) then + sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb + else + sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb + endif + else + sscalgrad2=0.0d0 + endif + return + end +c---------------------------------------------------------------------------- subroutine e_saxs(Esaxs_constr) implicit none include 'DIMENSIONS' @@ -11084,6 +11188,7 @@ c & sigma2CACA,sigma2CASC,sigma2SCCA,sigma2SCSC,expCACA,expCASC, & expSCCA,expSCSC,CASCgrad,SCCAgrad,SCSCgrad,aux,auxC,auxC1, & auxX,auxX1,CACAgrad,Cnorm + double precision sss2,ssgrad2,rrr,sscalgrad2,sscale2 double precision dist external dist c SAXS restraint penalty function @@ -11177,71 +11282,96 @@ c SC SC enddo ! k #else dijCACA=dist(i,j) - sigma2CACA=0.25d0/(restok(itype(j))**2+restok(itype(i))**2) + sigma2CACA=scal_rad**2*0.25d0/ + & (restok(itype(j))**2+restok(itype(i))**2) + + IF (saxs_cutoff.eq.0) THEN do k=1,nsaxs dk = distsaxs(k) expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2) Pcalc(k) = Pcalc(k)+expCACA + CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + do l=1,3 + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + enddo ! l + enddo ! k + ELSE + rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA) + do k=1,nsaxs + dk = distsaxs(k) +c write (2,*) "ijk",i,j,k + sss2 = sscale2(dijCACA,rrr,dk,0.3d0) + if (sss2.eq.0.0d0) cycle + ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0) + if (energy_dec) write(iout,'(a4,3i5,5f10.4)') + & 'saxs',i,j,k,dijCACA,rrr,dk,sss2,ssgrad2 + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2 + Pcalc(k) = Pcalc(k)+expCACA #ifdef DEBUG write(iout,*) "i j k Pcalc",i,j,Pcalc(k) #endif - CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + CACAgrad = -sigma2CACA*(dijCACA-dk)*expCACA+ + & ssgrad2*expCACA/sss2 do l=1,3 c CA CA aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA - PgradC(k,l,i) = PgradC(k,l,i)-aux - PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradC(k,l,i) = PgradC(k,l,i)+aux + PgradC(k,l,j) = PgradC(k,l,j)-aux enddo ! l enddo ! k + ENDIF #endif enddo ! j enddo ! iint enddo ! i #ifdef MPI if (nfgtasks.gt.1) then - call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION, - & MPI_SUM,king,FG_COMM,IERR) - if (fg_rank.eq.king) then + call MPI_AllReduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION, + & MPI_SUM,FG_COMM,IERR) +c if (fg_rank.eq.king) then do k=1,nsaxs Pcalc(k) = Pcalc_(k) enddo - endif - call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres, - & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) - if (fg_rank.eq.king) then - do i=1,nres - do l=1,3 - do k=1,nsaxs - PgradC(k,l,i) = PgradC_(k,l,i) - enddo - enddo - enddo - endif +c endif +c call MPI_AllReduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres, +c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR) +c if (fg_rank.eq.king) then +c do i=1,nres +c do l=1,3 +c do k=1,nsaxs +c PgradC(k,l,i) = PgradC_(k,l,i) +c enddo +c enddo +c enddo +c endif #ifdef ALLSAXS - call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres, - & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) - if (fg_rank.eq.king) then - do i=1,nres - do l=1,3 - do k=1,nsaxs - PgradX(k,l,i) = PgradX_(k,l,i) - enddo - enddo - enddo - endif +c call MPI_AllReduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres, +c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR) +c if (fg_rank.eq.king) then +c do i=1,nres +c do l=1,3 +c do k=1,nsaxs +c PgradX(k,l,i) = PgradX_(k,l,i) +c enddo +c enddo +c enddo +c endif #endif endif #endif -#ifdef MPI - if (fg_rank.eq.king) then -#endif Cnorm = 0.0d0 do k=1,nsaxs Cnorm = Cnorm + Pcalc(k) enddo - Esaxs_constr = dlog(Cnorm) +#ifdef MPI + if (fg_rank.eq.king) then +#endif + Esaxs_constr = dlog(Cnorm)-wsaxs0 do k=1,nsaxs - Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) + if (Pcalc(k).gt.0.0d0) + & Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) #ifdef DEBUG write (iout,*) "k",k," Esaxs_constr",Esaxs_constr #endif @@ -11249,6 +11379,11 @@ c CA CA #ifdef DEBUG write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr #endif +#ifdef MPI + endif +#endif + gsaxsC=0.0d0 + gsaxsX=0.0d0 do i=nnt,nct do l=1,3 auxC=0.0d0 @@ -11256,7 +11391,8 @@ c CA CA auxX=0.0d0 auxX1=0.d0 do k=1,nsaxs - auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) + if (Pcalc(k).gt.0) + & auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) auxC1 = auxC1+PgradC(k,l,i) #ifdef ALLSAXS auxX = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k) @@ -11272,7 +11408,7 @@ c * " gradX",wsaxs*(auxX - auxX1/Cnorm) enddo enddo #ifdef MPI - endif +c endif #endif return end