Add initial support for CASTEP file format.

src/recpot.c

@@ -273,9 +273,9 @@ static int _to_vr(pspio_potential_t *pot, const double *qs, const double *vq, si
                                   PSPIO_MESH_LOG2, LOG_A, LOG_B);
   rs = pspio_mesh_get_r(pot->v->mesh);
 
-  pot->v->f[0] = !pot->qn->l ? _trapz(vq, qs, nq) : 0.;
+  pot->v->f[0] = !pot->qn->l ? 0.5 * _trapz(vq, qs, nq) : 0.;
   for (int ir = 1; ir < pot->v->mesh->np; ir ++)
-    pot->v->f[ir] = _trapz_jl(pot->qn->l, vq, qs, rs[ir], nq);
+    pot->v->f[ir] = 0.5 * _trapz_jl(pot->qn->l, vq, qs, rs[ir], nq);
   SUCCEED_OR_RETURN( pspio_interp_init(pot->v->f_interp,
                                        pot->v->mesh, pot->v->f) );
   for (int ir = 1; ir < pot->v->mesh->np; ir ++)
@@ -314,9 +314,9 @@ static int _to_projs(pspio_projector_t *p, const double *qs, const double *vq, s
                                   PSPIO_MESH_LOG2, LOG_A, LOG_B);
   rs = pspio_mesh_get_r(p->proj->mesh);
 
-  p->proj->f[0] = !p->qn->l ? _trapz(vq, qs, nq) : 0.;
+  p->proj->f[0] = !p->qn->l ? 0.5 * _trapz(vq, qs, nq) : 0.;
   for (int ir = 1; ir < p->proj->mesh->np; ir ++)
-    p->proj->f[ir] = _trapz_jl(p->qn->l, vq, qs, rs[ir], nq);
+    p->proj->f[ir] = 0.5 * _trapz_jl(p->qn->l, vq, qs, rs[ir], nq);
   SUCCEED_OR_RETURN( pspio_interp_init(p->proj->f_interp,
                                        p->proj->mesh, p->proj->f) );
   for (int ir = 1; ir < p->proj->mesh->np; ir ++)
@@ -366,8 +366,8 @@ int pspio_recpot_read(FILE *fp, pspio_pspdata_t *pspdata)
   SUCCEED_OR_RETURN( recpot_to_projs(&recpot, pspdata) );
 
   z = round(pspdata->vlocal->v->f[pspdata->vlocal->v->mesh->np - 1] * pspdata->vlocal->v->mesh->r[pspdata->vlocal->v->mesh->np - 1]);
-  pspio_pspdata_set_zvalence(pspdata, - z * 0.5);
-  pspio_pspdata_set_nelvalence(pspdata, - z * 0.5);
+  pspio_pspdata_set_zvalence(pspdata, - z);
+  pspio_pspdata_set_nelvalence(pspdata, - z);
 
   return PSPIO_SUCCESS;
 }
@@ -404,9 +404,9 @@ int pspio_recpot_write(FILE *fp, const pspio_pspdata_t *pspdata)
   snprintf(line, PSPIO_STRLEN_LINE, "%18.8f\n", dq * nq);
   _write_str(fp, line);
 
-  vlq[0] = _trapz_r(pspdata->vlocal->v->f, pspdata->vlocal->v->mesh->r, 2., pspdata->vlocal->v->mesh->np);
+  vlq[0] = 2. * _trapz_r(pspdata->vlocal->v->f, pspdata->vlocal->v->mesh->r, 2., pspdata->vlocal->v->mesh->np);
   for (int iq = 1; iq < nq + 1; iq ++)
-    vlq[iq] = _trapz_qr(pspdata->vlocal->v->f, pspdata->vlocal->v->mesh->r, dq * iq * bohrad, 20., pspdata->vlocal->v->mesh->np);
+    vlq[iq] = 2. * _trapz_qr(pspdata->vlocal->v->f, pspdata->vlocal->v->mesh->r, dq * iq * bohrad, 20., pspdata->vlocal->v->mesh->np);
   for (int iq = 0; iq < nq + 1; iq += 3) {
     snprintf(line, PSPIO_STRLEN_LINE, " %19.10f %19.10f %19.10f\n", vlq[iq], vlq[iq + 1], vlq[iq + 2]);
     _write_str(fp, line);
@@ -418,9 +418,9 @@ int pspio_recpot_write(FILE *fp, const pspio_pspdata_t *pspdata)
     _write_str(fp, line);
     snprintf(line, PSPIO_STRLEN_LINE, "%20.10f\n", pspio_projector_get_energy(proj) * ryd);
     _write_str(fp, line);
-    vlq[0] = !pspio_qn_get_l(proj->qn) ? _trapz_q(proj->proj->f, proj->proj->mesh->r, 20., proj->proj->mesh->np): 0.;
+    vlq[0] = !pspio_qn_get_l(proj->qn) ? 2. * _trapz_q(proj->proj->f, proj->proj->mesh->r, 20., proj->proj->mesh->np): 0.;
     for (int iq = 1; iq < nq + 1; iq ++)
-      vlq[iq] = _trapz_jlq(pspio_qn_get_l(proj->qn), proj->proj->f, proj->proj->mesh->r, dq * iq * bohrad, 2., proj->proj->mesh->np);
+      vlq[iq] = 2. * _trapz_jlq(pspio_qn_get_l(proj->qn), proj->proj->f, proj->proj->mesh->r, dq * iq * bohrad, 2., proj->proj->mesh->np);
     for (int iq = 0; iq < nq + 1; iq += 3) {
       snprintf(line, PSPIO_STRLEN_LINE, " %19.10f %19.10f %19.10f\n", vlq[iq], vlq[iq + 1], vlq[iq + 2]);
       _write_str(fp, line);