anchor.diff:          call nwpw_GVT4(aa,bb,cc,dd,ee,ff,alpha,alpha,
        >                 xrud,zrud,gammaud,GGud,dGud_dxud,dGud_dzud)
   
  -        eudc        = GGud*eud0c
  -        dGud_dnup   = dGud_dxud*dxup_dnup + dGud_dzud*dzup_dnup
  -        dGud_dagrup = dGud_dxud*dxup_dagrup 
  +        dGud_dnup   = dGud_dxud*dxup_dnup   + dGud_dzud*dzup_dnup
  +        dGud_dndn   = dGud_dxud*dxdn_dndn   + dGud_dzud*dzdn_dndn
  +        dGud_dagrup = dGud_dxud*dxup_dagrup
  +        dGud_dagrdn = dGud_dxud*dxdn_dagrdn
           dGud_dtauup = dGud_dzud*dzup_dtauup
  -        dGud_dndn   = dGud_dxud*dxdn_dndn + dGud_dzud*dzdn_dndn
  -        dGud_dagrdn = dGud_dxud*dxdn_dagrdn 
           dGud_dtaudn = dGud_dzud*dzdn_dtaudn
   
  -        !fudc = n*GGud*eud0c
  -        fud_nupc    = GGud*eud0c 
  -     &                + nup*(dGud_dnup*eud0c + GGud*fud0_nupc)
  -        fud_dnupc   = nup*eud0c*dGud_dagrup
  -        fud_dtauupc = nup*eud0c*dGud_dtauup
  -        fud_ndnc    = GGud*eud0c 
  -     &                + ndn*(dGud_dndn*eud0c + GGud*fud0_ndnc)
  -        fud_dndnc   = ndn*eud0c*dGud_dagrdn
  -        fud_dtaudnc = ndn*eud0c*dGud_dtaudn
  -
  -        ec       = euuc + eddc + eudc
  -        fnupc    = fuu_nupc + fud_nupc
  -        fdnupc   = fuu_dnupc + fud_dnupc
  -        fdtauupc = fuu_dtauupc + fud_dtauupc 
  -        fndnc    = fdd_ndnc + fud_ndnc
  -        fdndnc   = fdd_dndnc + fud_dndnc
  -        fdtaudnc = fdd_dtaudnc + fud_dtaudnc 
  -        
  +        eudc   = GGud*eud1c
  +
  +        dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn 
  +     &              + deud0c_dzeta*dzeta_dnup) - euu0c - nup*deuu0c_dnup
  +        dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn 
  +     &              + deud0c_dzeta*dzeta_dndn) - edd0c - ndn*dedd0c_dndn
  +
  +        dfudc_dnup   = fud1c*dGud_dnup    + GGud*dfud1c_dnup
  +        dfudc_dndn   = fud1c*dGud_dndn    + GGud*dfud1c_dndn
  +        dfudc_dagrup = fud1c*dGud_dagrup
  +        dfudc_dagrdn = fud1c*dGud_dagrdn
  +        dfudc_dtauup = fud1c*dGud_dtauup
  +        dfudc_dtaudn = fud1c*dGud_dtaudn
  +
  +        ec       = eudc          + (euuc*nup  + eddc*ndn)/n 
  +        fnupc    = dfudc_dnup    + dfuuc_dnup
  +        fndnc    = dfudc_dndn    + dfddc_dndn
  +        fdnupc   = dfudc_dagrup  + dfuuc_dagrup
  +        fdndnc   = dfudc_dagrdn  + dfddc_dagrdn
  +        fdtauupc = dfudc_dtauup  + dfuuc_dtauup
  +        fdtaudnc = dfudc_dtaudn  + dfddc_dtaudn
  +
           xce(i)     = x_parameter*ex       + c_parameter*ec
           fn(i,1)    = x_parameter*fnupx    + c_parameter*fnupc
           fn(i,2)    = x_parameter*fndnx    + c_parameter*fndnc
  
positive.code:         ee    = -1.282732d-4
          ff    =  3.574822d-4
          alpha =  0.00186726d0
  
  ********* Spin-up *********
          gammaup = 1.0d0 + alpha*(xup + zup)
          xrup    = xup/gammaup
          zrup    = zup/gammaup
  
          call nwpw_GVT4(aa,bb,cc,dd,ee,ff,alpha,alpha,
       >                 xrup,zrup,gammaup,GGup,dGup_dxup,dGup_dzup)
  
          eupx     = nup_13*GGup
          fnupx    = frthrd*nup_13*GGup 
       &           + nup_43*(dGup_dxup*dxup_dnup + dGup_dzup*dzup_dnup)
          fdnupx   = nup_43*(dGup_dxup*dxup_dagrup)
          fdtauupx = nup_43*(dGup_dzup*dzup_dtauup) 
  ********* Spin-down *********
          gammadn = 1.0d0 + alpha*(xdn + zdn)
          xrdn    = xdn/gammadn
          zrdn    = zdn/gammadn
  
          call nwpw_GVT4(aa,bb,cc,dd,ee,ff,alpha,alpha,
       >                 xrdn,zrdn,gammadn,GGdn,dGdn_dxdn,dGdn_dzdn)
  
          ednx     = nup_13*GGdn
          fndnx    = frthrd*ndn_13*GGdn 
       &           + ndn_43*(dGdn_dxdn*dxdn_dndn + dGdn_dzdn*dzdn_dndn)
          fdndnx   = ndn_43*(dGdn_dxdn*dxdn_dagrdn)
          fdtaudnx = ndn_43*(dGdn_dzdn*dzdn_dtaudn) 
   
          ex = eupx + ednx
  ********* VS98 Correlation *********
  ********* Same-spin *********
          aa    =  3.270912d-1
          bb    = -3.228915d-2
          cc    = -2.942406d-2
  
negative.code:         ee    = -5.451559d-3
          ff    =  1.577575d-2
          alpha =  0.00515088d0
  *       ***** alpha-alpha *****
          gammaup = 1.0d0 + alpha*(xup + zup)
          xrup    = xup/gammaup
          zrup    = zup/gammaup
  
          call nwpw_vs98c_ss(aa,bb,cc,dd,ee,ff,alpha,nup,agrup,tauup,
       >               xup,zup,gammaup,xrup,zrup,
       >               dxup_dnup,dxup_dagrup,dzup_dnup,dzup_dtauup,
       >               euu0c,deuu0c_dnup,
       >               euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup)
  
  *       ***** beta-beta *****
          gammadn = 1.0d0 + alpha*(xdn + zdn)
          xrdn    = xdn/gammadn
          zrdn    = zdn/gammadn
  
          call nwpw_vs98c_ss(aa,bb,cc,dd,ee,ff,alpha,ndn,agrdn,taudn,
       >               xdn,zdn,gammadn,xrdn,zrdn,
       >               dxdn_dndn,dxdn_dagrdn,dzdn_dndn,dzdn_dtaudn,
       >               edd0c,dedd0c_dndn,
       >               eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn)
  
  *       ***** opposite-spin *****
          aa    =  7.035010d-1
          bb    =  7.694574d-3
          cc    =  5.152765d-2
  