Changeset 600

Timestamp:

04/16/11 19:38:57 (2 years ago)

Author:

micael

Message:

• 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:

• trunk/src/ps_io.F90 (modified) (3 diffs)

trunk/src/ps_io.F90

@@ -687 +687 @@
     !Mesh
     call mesh_null(new_m)
-    call mesh_init(new_m, MESH_LOG1, MESH_CUBIC_SPLINE, info%m%r(1), info%m%r(info%m%np), np=info%m%np)
+    call mesh_init(new_m, MESH_LOG1, MESH_FINITE_DIFF, info%m%r(1), info%m%r(info%m%np), np=info%m%np, fd_order=4)
 
     !Header
@@ -758 +758 @@
     integer, intent(in) :: unit
 
-    integer :: values(8)
+    integer :: i, values(8)
+    real(R8) :: rc
     character(60) :: header
 
@@ -782 +783 @@
           ixc_to_abinit(info%ixc), maxval(info%psp_l), info%kb_l_local, info%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*mesh_integrate(info%m, info%rho_core), 'rchrg,fchrg,qchrg'
     else
       write(unit,'(3F20.14,T64,A)') 0.0, 0.0, 0.0, 'rchrg,fchrg,qchrg'