Changeset 604

Timestamp:

06/26/11 18:03:07 (2 years ago)

Author:

micael

Message:

Backport of #600

• Bug fixed: the definition of the core charge radius in the Abinit file format is in fact the radius at which the core density vanishes.

File:

• branches/1.1.x/src/ps_io.F90 (modified) (2 diffs)

branches/1.1.x/src/ps_io.F90

@@ -718 +718 @@
     integer, intent(in) :: unit
 
-    integer :: values(8)
+    integer :: i, values(8)
+    real(R8) :: rc
     type(mesh_type) :: new_m
     character(60) :: header
@@ -745 +746 @@
           ixc_to_abinit(info%ixc), maxval(info%psp_l), info%kb_l_local, new_m%np, 0.0
     if (info%nlcc) then
-      write(unit,'(3F20.14,T64,A)') info%nlcc_rc, 1.0, 0.0, 'rchrg,fchrg,qchrg'
+      do i = info%m%np, 1, -1 
+        if (info%rho_core(i) > 1.e-6) then 
+          rc = info%m%r(i) 
+          exit 
+        end if
+      end do
+      write(unit,'(3F20.14,T64,A)') rc, 1.0, &
+           M_FOUR*M_PI*spline_integ(info%m%np, info%m%r, info%rho_core*info%m%r**2, M_ZERO, info%m%r(info%m%np), 3), &
+           'rchrg,fchrg,qchrg'
     else
       write(unit,'(3F20.14,T64,A)') 0.0, 0.0, 0.0, 'rchrg,fchrg,qchrg'