LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 12:23:48 +01:00
Code Issues Releases Wiki Activity

Revert input coordinates

Browse Source
This commit is contained in:
Anthony Scemama 2017-12-19 11:49:48 +01:00
parent ac37a499d3
commit 17e0518410
6 changed files with 235 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
plugins/Symmetry/Symmetry.main.irp.f
View File

@ -6,13 +6,13 @@ program Symmetry
integer :: i, j
character*8 :: sym
print *, 'Molecule is linear:', molecule_is_linear
print *, 'Has center of inversion:', molecule_has_center_of_inversion
print *, 'Has S2n improper rotation:', molecule_has_improper_rotation
print *, 'Symmetry rotation axis:', sym_rotation_axis(:)
print *, 'Group:'//point_group
print *, 'Symmetry irreps', sym_irrep(1:n_irrep)
print *, 'Symmetry operations', sym_operation(1:n_irrep)
print *, 'Molecule is linear: ', molecule_is_linear
print *, 'Has center of inversion: ', molecule_has_center_of_inversion
print *, 'Has S2n improper rotation: ', molecule_has_improper_rotation
print *, 'Symmetry rotation axis: ', sym_rotation_axis(:)
print *, 'Group: '//point_group
print *, 'Symmetry irreps : ', sym_irrep(1:n_irrep)
print *, 'Symmetry operations : ', sym_operation(1:n_irrep)
print *, 'Character table'
do i=1,n_irrep
print *, character_table(i,:)

38
plugins/Symmetry/aos.irp.f
View File

@ -7,8 +7,8 @@ BEGIN_PROVIDER [ double precision, sym_box, (3,2) ]
sym_box(:,:) = 0.d0
do xyz=1,3
do i=1,nucl_num
sym_box(xyz,1) = min(sym_box(xyz,1), nucl_coord(i,xyz))
sym_box(xyz,2) = max(sym_box(xyz,2), nucl_coord(i,xyz))
sym_box(xyz,1) = min(sym_box(xyz,1), nucl_coord_sym(i,xyz))
sym_box(xyz,2) = max(sym_box(xyz,2), nucl_coord_sym(i,xyz))
enddo
enddo
sym_box(:,1) = sym_box(:,1) - 2.d0
@ -42,24 +42,28 @@ subroutine compute_sym_ao_values(sym_points, n_sym_points, result)
double precision, intent(out) :: result(n_sym_points, ao_num)
integer :: i, j
double precision :: x, y, z
double precision :: x2, y2, z2
result (:,:) = 0.d0
do j=1,ao_num
do i=1,n_sym_points
x = sym_points(1,i) - nucl_coord_transp(1,ao_nucl(j))
y = sym_points(2,i) - nucl_coord_transp(2,ao_nucl(j))
z = sym_points(3,i) - nucl_coord_transp(3,ao_nucl(j))
x = x**ao_power(j,1)
y = y**ao_power(j,2)
z = z**ao_power(j,3)
! result(i,j) = x*y*z*exp(-(x*x+y*y+z*z))
result(i,j) = x*y*z
if (result(i,j) > 0.d0) then
result(i,j) = 1.d0
else if (result(i,j) < 0.d0) then
result(i,j) = -1.d0
else
result(i,j) = 0.d0
endif
x = sym_points(1,i) - nucl_coord_sym_transp(1,ao_nucl(j))
y = sym_points(2,i) - nucl_coord_sym_transp(2,ao_nucl(j))
z = sym_points(3,i) - nucl_coord_sym_transp(3,ao_nucl(j))
x2 = x*sym_molecule_rotation_inv(1,1) + y*sym_molecule_rotation_inv(2,1) + z*sym_molecule_rotation_inv(3,1)
y2 = x*sym_molecule_rotation_inv(1,2) + y*sym_molecule_rotation_inv(2,2) + z*sym_molecule_rotation_inv(3,2)
z2 = x*sym_molecule_rotation_inv(1,3) + y*sym_molecule_rotation_inv(2,3) + z*sym_molecule_rotation_inv(3,3)
x = x2**ao_power(j,1)
y = y2**ao_power(j,2)
z = z2**ao_power(j,3)
result(i,j) = x*y*z*exp(-(x*x+y*y+z*z))
! result(i,j) = x*y*z
! if (result(i,j) > 0.d0) then
! result(i,j) = 1.d0
! else if (result(i,j) < 0.d0) then
! result(i,j) = -1.d0
! else
! result(i,j) = 0.d0
! endif
enddo
enddo

21
plugins/Symmetry/find_sym.irp.f
View File

@ -20,7 +20,7 @@ BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ]
found = .False.
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j)
point(:) = nucl_coord_sym_transp(:,i) + nucl_coord_sym_transp(:,j)
if (u_dot_u(point,3) < 1.d-5) then
found = .True.
exit
@ -52,10 +52,10 @@ BEGIN_PROVIDER [ integer, sym_rotation_axis, (3) ]
sym_rotation_axis(iaxis) = iorder
do i=1,nucl_num
found = .False.
call sym_apply_rotation(dble(iorder),iaxis,nucl_coord_transp(1,i),point)
call sym_apply_rotation(dble(iorder),iaxis,nucl_coord_sym_transp(1,i),point)
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point2(:) = nucl_coord_transp(:,j) - point(:)
point2(:) = nucl_coord_sym_transp(:,j) - point(:)
if (u_dot_u(point2,3) < 1.d-5) then
found = .True.
exit
@ -148,10 +148,10 @@ BEGIN_PROVIDER [ logical, molecule_has_improper_rotation ]
molecule_has_improper_rotation = .True.
do i=1,nucl_num
found = .False.
call sym_apply_improper_rotation(dble(iorder),iaxis,nucl_coord_transp(1,i),point)
call sym_apply_improper_rotation(dble(iorder),iaxis,nucl_coord_sym_transp(1,i),point)
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point2(:) = nucl_coord_transp(:,j) - point(:)
point2(:) = nucl_coord_sym_transp(:,j) - point(:)
if (u_dot_u(point2,3) < 1.d-5) then
found = .True.
exit
@ -179,7 +179,7 @@ BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ]
found = .False.
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j)
point(:) = nucl_coord_sym_transp(:,i) + nucl_coord_sym_transp(:,j)
if (u_dot_u(point,3) < 1.d-5) then
found = .True.
exit
@ -208,11 +208,11 @@ BEGIN_PROVIDER [ logical, molecule_has_sigma_plane, (3) ]
molecule_has_sigma_plane(iaxis) = .True.
do i=1,nucl_num
found = .False.
point(:) = nucl_coord_transp(:,i)
point(:) = nucl_coord_sym_transp(:,i)
point(iaxis) = -point(iaxis)
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point2(:) = nucl_coord_transp(:,j) - point(:)
point2(:) = nucl_coord_sym_transp(:,j) - point(:)
if (u_dot_u(point2,3) < 1.d-5) then
found = .True.
exit
@ -233,7 +233,7 @@ BEGIN_PROVIDER [ character*16, point_group ]
! Point group of the molecule
END_DOC
character*2, save :: i_to_a(24) = (/ ' 1', ' 2', ' 3', ' 4', ' 5', ' 6', ' 7', ' 8', ' 9', &
character*2, save :: i_to_a(24) = (/ '1 ', '2 ', '3 ', '4 ', '5 ', '6 ', '7 ', '8 ', '9 ', &
'10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '20', &
'21', '22', '23', '24' /)
point_group = 'C1'
@ -366,7 +366,7 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ]
integer :: iangle, n_sym_points
double precision :: angle
integer :: iop, imo, ipoint, l, i
double precision :: sym_operations_on_mos(n_irrep)
double precision :: sym_operations_on_mos(mo_tot_num)
logical :: possible_irrep(n_irrep,mo_tot_num)
n_sym_points = 10000
@ -443,6 +443,7 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ]
sym_operations_on_mos(imo) += x
enddo
sym_operations_on_mos(imo) *= 1.d0/n_sym_points
print *, iop, imo, sym_operations_on_mos(imo)
if (dabs(sym_operations_on_mos(imo)-1.d0) < 1.d-2) then
sym_operations_on_mos(imo)=1.d0
else if (dabs(sym_operations_on_mos(imo)+1.d0) < 1.d-2) then

92
plugins/Symmetry/nuclei.irp.f Normal file
View File

@ -0,0 +1,92 @@
BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates in standard orientation
END_DOC
if (mpi_master) then
integer :: i
do i=1,nucl_num
nucl_coord_sym(i,1) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + &
(nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + &
(nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1)
nucl_coord_sym(i,2) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + &
(nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + &
(nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2)
nucl_coord_sym(i,3) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + &
(nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + &
(nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3)
enddo
character*(64), parameter :: f = '(A16, 4(1X,F12.6))'
character*(64), parameter :: ft= '(A16, 4(1X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(output_Nuclei)
write(output_Nuclei,'(A)') ''
write(output_Nuclei,'(A)') 'Nuclear Coordinates in standard orientation (Angstroms)'
write(output_Nuclei,'(A)') '======================================================='
write(output_Nuclei,'(A)') ''
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,*) &
' Atom Charge X Y Z '
write(output_Nuclei,ft) &
'================','============','============','============','============'
do i=1,nucl_num
write(output_Nuclei,f) nucl_label(i), nucl_charge(i), &
nucl_coord_sym(i,1)*a0, &
nucl_coord_sym(i,2)*a0, &
nucl_coord_sym(i,3)*a0
enddo
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,'(A)') ''
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( nucl_coord_sym, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read nucl_coord_sym with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord_sym_transp, (3,nucl_num) ]
implicit none
BEGIN_DOC
! Transposed array of nucl_coord
END_DOC
integer :: i, k
nucl_coord_sym_transp = 0.d0
do i=1,nucl_num
nucl_coord_sym_transp(1,i) = nucl_coord_sym(i,1)
nucl_coord_sym_transp(2,i) = nucl_coord_sym(i,2)
nucl_coord_sym_transp(3,i) = nucl_coord_sym(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, sym_molecule_rotation, (3,3) ]
implicit none
BEGIN_DOC
! Rotation of the molecule to go from input orientation to standard orientation
END_DOC
call find_rotation(nucl_coord, size(nucl_coord,1), nucl_coord_sym, 3, sym_molecule_rotation, 3)
END_PROVIDER
BEGIN_PROVIDER [ double precision, sym_molecule_rotation_inv, (3,3) ]
implicit none
BEGIN_DOC
! Rotation of the molecule to go from standard orientation to input orientation
END_DOC
call find_rotation(nucl_coord_sym, size(nucl_coord_sym,1), nucl_coord, 3, sym_molecule_rotation_inv, 3)
END_PROVIDER

12
src/Nuclei/inertia.irp.f
View File

@ -6,12 +6,12 @@ BEGIN_PROVIDER [ double precision, inertia_tensor, (3,3) ]
integer :: i,j,k
inertia_tensor = 0.d0
do k=1,nucl_num
inertia_tensor(1,1) += element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2)
inertia_tensor(2,2) += element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2)
inertia_tensor(3,3) += element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,2)-center_of_mass(2))**2)
inertia_tensor(1,2) -= element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1)) * (nucl_coord_input(k,2)-center_of_mass(2)) )
inertia_tensor(1,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1)) * (nucl_coord_input(k,3)-center_of_mass(3)) )
inertia_tensor(2,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2)) * (nucl_coord_input(k,3)-center_of_mass(3)) )
inertia_tensor(1,1) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,2)-center_of_mass(2))**2 + (nucl_coord(k,3)-center_of_mass(3))**2)
inertia_tensor(2,2) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1))**2 + (nucl_coord(k,3)-center_of_mass(3))**2)
inertia_tensor(3,3) += element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1))**2 + (nucl_coord(k,2)-center_of_mass(2))**2)
inertia_tensor(1,2) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1)) * (nucl_coord(k,2)-center_of_mass(2)) )
inertia_tensor(1,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,1)-center_of_mass(1)) * (nucl_coord(k,3)-center_of_mass(3)) )
inertia_tensor(2,3) -= element_mass(int(nucl_charge(k))) * ((nucl_coord(k,2)-center_of_mass(2)) * (nucl_coord(k,3)-center_of_mass(3)) )
enddo
inertia_tensor(2,1) = inertia_tensor(1,2)
inertia_tensor(3,1) = inertia_tensor(1,3)

223
src/Nuclei/nuclei.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ]
BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ]
implicit none
BEGIN_DOC
@ -8,7 +8,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ]
if (mpi_master) then
double precision, allocatable :: buffer(:,:)
nucl_coord_input = 0.d0
nucl_coord = 0.d0
allocate (buffer(nucl_num,3))
buffer = 0.d0
logical :: has
@ -22,7 +22,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ]
do i=1,3
do j=1,nucl_num
nucl_coord_input(j,i) = buffer(j,i)
nucl_coord(j,i) = buffer(j,i)
enddo
enddo
deallocate(buffer)
@ -31,65 +31,6 @@ BEGIN_PROVIDER [ double precision, nucl_coord_input, (nucl_num,3) ]
character*(64), parameter :: ft= '(A16, 4(1X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(output_Nuclei)
write(output_Nuclei,'(A)') ''
write(output_Nuclei,'(A)') 'Input Nuclear Coordinates (Angstroms)'
write(output_Nuclei,'(A)') '====================================='
write(output_Nuclei,'(A)') ''
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,*) &
' Atom Charge X Y Z '
write(output_Nuclei,ft) &
'================','============','============','============','============'
do i=1,nucl_num
write(output_Nuclei,f) nucl_label(i), nucl_charge(i), &
nucl_coord_input(i,1)*a0, &
nucl_coord_input(i,2)*a0, &
nucl_coord_input(i,3)*a0
enddo
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,'(A)') ''
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( nucl_coord_input, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read nucl_coord_input with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates in standard orientation
END_DOC
if (mpi_master) then
integer :: i
do i=1,nucl_num
nucl_coord(i,1) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + &
(nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + &
(nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1)
nucl_coord(i,2) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + &
(nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + &
(nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2)
nucl_coord(i,3) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + &
(nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + &
(nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3)
enddo
character*(64), parameter :: f = '(A16, 4(1X,F12.6))'
character*(64), parameter :: ft= '(A16, 4(1X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(output_Nuclei)
write(output_Nuclei,'(A)') ''
write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)'
@ -205,35 +146,36 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
END_DOC
PROVIDE mpi_master nucl_coord nucl_charge nucl_num
if (disk_access_nuclear_repulsion.EQ.'Read') then
logical :: has
if (mpi_master) then
call ezfio_has_nuclei_nuclear_repulsion(has)
if (has) then
call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion)
else
print *, 'nuclei/nuclear_repulsion not found in EZFIO file'
stop 1
if (mpi_master) then
if (disk_access_nuclear_repulsion.EQ.'Read') then
logical :: has
if (mpi_master) then
call ezfio_has_nuclei_nuclear_repulsion(has)
if (has) then
call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion)
else
print *, 'nuclei/nuclear_repulsion not found in EZFIO file'
stop 1
endif
print*, 'Read nuclear_repulsion'
endif
print*, 'Read nuclear_repulsion'
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read nuclear_repulsion with MPI'
endif
IRP_ENDIF
else
integer :: k,l
double precision :: Z12, r2, x(3)
nuclear_repulsion = 0.d0
do l = 1, nucl_num
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read nuclear_repulsion with MPI'
endif
IRP_ENDIF
else
integer :: k,l
double precision :: Z12, r2, x(3)
nuclear_repulsion = 0.d0
do l = 1, nucl_num
do k = 1, nucl_num
if(k == l) then
cycle
@ -245,45 +187,76 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3)
nuclear_repulsion += Z12/dsqrt(r2)
enddo
enddo
nuclear_repulsion *= 0.5d0
end if
call write_time(output_Nuclei)
call write_double(output_Nuclei,nuclear_repulsion, &
'Nuclear repulsion energy')
if (disk_access_nuclear_repulsion.EQ.'Write') then
if (mpi_master) then
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
enddo
nuclear_repulsion *= 0.5d0
end if
call write_time(output_Nuclei)
call write_double(output_Nuclei,nuclear_repulsion, &
'Nuclear repulsion energy')
if (disk_access_nuclear_repulsion.EQ.'Write') then
if (mpi_master) then
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
endif
endif
endif
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( element_name, size(element_name)*4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read element_name with MPI'
endif
call MPI_BCAST( element_mass, size(element_mass), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read element_name with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ character*(4), element_name, (0:128)]
BEGIN_PROVIDER [ character*(4), element_name, (0:128)]
&BEGIN_PROVIDER [ double precision, element_mass, (0:128) ]
BEGIN_DOC
! Array of the name of element, sorted by nuclear charge (integer)
END_DOC
integer :: iunit
integer, external :: getUnitAndOpen
character*(128) :: filename
call getenv('QP_ROOT',filename)
filename = trim(filename)//'/data/list_element.txt'
iunit = getUnitAndOpen(filename,'r')
element_mass(:) = 0.d0
do i=0,128
write(element_name(i),'(I4)') i
enddo
character*(80) :: buffer, dummy
do
read(iunit,'(A80)',end=10) buffer
read(buffer,*) i ! First read i
read(buffer,*) i, element_name(i), dummy, element_mass(i)
enddo
10 continue
close(10)
BEGIN_DOC
! Array of the name of element, sorted by nuclear charge (integer)
END_DOC
integer :: iunit
integer, external :: getUnitAndOpen
character*(128) :: filename
if (mpi_master) then
call getenv('QP_ROOT',filename)
filename = trim(filename)//'/data/list_element.txt'
iunit = getUnitAndOpen(filename,'r')
element_mass(:) = 0.d0
do i=0,128
write(element_name(i),'(I4)') i
enddo
character*(80) :: buffer, dummy
do
read(iunit,'(A80)',end=10) buffer
read(buffer,*) i ! First read i
read(buffer,*) i, element_name(i), dummy, element_mass(i)
enddo
10 continue
close(10)
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( element_name, size(element_name)*4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read element_name with MPI'
endif
call MPI_BCAST( element_mass, size(element_mass), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read element_name with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ double precision, center_of_mass, (3) ]
@ -297,7 +270,7 @@ BEGIN_PROVIDER [ double precision, center_of_mass, (3) ]
s = 0.d0
do i=1,nucl_num
do j=1,3
center_of_mass(j) += nucl_coord_input(i,j)* element_mass(int(nucl_charge(i)))
center_of_mass(j) += nucl_coord(i,j)* element_mass(int(nucl_charge(i)))
enddo
s += element_mass(int(nucl_charge(i)))
enddo