diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f index 78bf126e..1d438af5 100644 --- a/plugins/Symmetry/Symmetry.main.irp.f +++ b/plugins/Symmetry/Symmetry.main.irp.f @@ -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,:) diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f index e4352948..e0aade1e 100644 --- a/plugins/Symmetry/aos.irp.f +++ b/plugins/Symmetry/aos.irp.f @@ -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 diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f index 80e2fccd..95bb4e6c 100644 --- a/plugins/Symmetry/find_sym.irp.f +++ b/plugins/Symmetry/find_sym.irp.f @@ -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 diff --git a/plugins/Symmetry/nuclei.irp.f b/plugins/Symmetry/nuclei.irp.f new file mode 100644 index 00000000..caa4781e --- /dev/null +++ b/plugins/Symmetry/nuclei.irp.f @@ -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 + diff --git a/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f index 97e32711..d2cb7cf8 100644 --- a/src/Nuclei/inertia.irp.f +++ b/src/Nuclei/inertia.irp.f @@ -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) diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 0d14ae7e..7e4de39f 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -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