nuclei_mod.f95
nuclei_mod.f95 ! -*- f90 -*-
!-----------------------------------------------------------------------------
! This file is part of QmcMol A GNU QMC program for molecules.
! Copyright (C) 2002 CNRS - U.P.M.C. Paris 6 - FRANCE
! R. Assaraf, M. Caffarel, F. Colonna, X. Krokidis, B. Levy,
! P. Pernod, and P. Reinhardt - qmcmol@lct.jussieu.fr
! http://www.lct.jussieu.fr/QmcMol
! QmcMol is free software; you can redistribute it and/or modify it
! under the terms of the GNU General Public License as published by the Free
! Software Foundation; either version 2, or (at your option) any later
! version. QmcMol is distributed in the hope that it will be useful,but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
! for more details. You should have received a copy of the GNU General
! Public License along with QmcMol; see the file COPYING. If not,
! write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
! Boston, MA 02111-1307, USA.
!-----------------------------------------------------------------------------
module nuclei_mod
use basic_types_tools_mod
use control_mod
use arrays_mod
use mendeleiev_mod
use parser_mod
use system_data_mod
! Revision: 13 Fev 2001: stemp, string => local
! Revision: 25 Sep 2002: max_len_atom_symbol moved to control_mod
! Revision: 21 Dec 2002: atoms_print_res
type atom
real(dp), pointer :: coordinates(:) => NULL ()
real(dp) :: charge, mass
character(len=max_len_pseudo_name) :: pseudo_name
character(len=max_len_atom_name) :: name
character(len=max_len_atom_symbol) :: symbol
end type atom
type (atom), pointer :: atoms(:) => NULL ()
! Improve NON pas de cur_atom
integer(i4b) :: cur_atom = -1
character(len=max_string_len) :: symmetry_type
real(dp) :: nuclear_potential_energy
real(dp), pointer :: nuclear_dipole_moment(:) => NULL()
real(dp), pointer :: center_of_mass(:) => NULL()
integer(i4b), pointer :: symmetry_array(:) => NULL()
contains
! =========================================================================
subroutine nuclei_menu
! -------------------------------------------------------------------------
!
! Name : nuclei_menu
!
! Description : input nuclei
!
! Authors:X. Krokidis, F. Colonna
! Date :26 Sep 2000
! -------------------------------------------------------------------------
implicit none
! local:
integer(i4b) :: memory_siz
character(len=max_string_len_routine_name) :: l_here
logical(bool) :: done
logical(bool) :: l_print
! begin:
l_here ='nuclei_menu'
call enter (l_here)
! initialize:
l_print = .false.
done = .false.
call next_command_get ('nuclei')
do while (.not.done)
call next_word_get ('nuclei')
if(trim(cur_word) .eq. 'help') then
write (logf,*) 'nuclei-h: nuclei definitions'
write (logf,*) 'nuclei-h: syntax:'
write (logf,*) 'nuclei-h: nuclei'
write (logf,*) 'nuclei-h: help'
write (logf,*) 'nuclei-h: ato{m} [command]'
write (logf,*) 'nuclei-h: pri{nt} = [boolean] results file'
write (logf,*) 'nuclei-h: sym{metry} [command]'
write (logf,*) 'nuclei-h: end'
write (logf,*) 'nuclei-h: example'
write (logf,*) 'nuclei-h: nuclei atom ... end print=atoms end_nuc'
else if(cur_word(1:3) .eq. 'ato') then
call atom_menu
else if(token_is('sym', list_of_integer)) then
call copy (symmetry_array, list_of_integer)
else if(token_is('pri', l_print)) then
else
call menu_exit (l_here, done)
end if
end do
message ='in routine '*l_here*' cur_atom='*cur_atom*CR* &
' size(atoms)='*size(atoms)
call ensure (message, cur_atom .eq. size(atoms))
if(debug) call atoms_print
! check
call atoms_check
if(l_print) then
call atoms_print
call atoms_print_res ('atoms information')
end if
memory_siz = (2+max_len_pseudo_name+max_len_atom_name &
+2*dp+dp*size(atoms(1)%coordinates))*size(atoms)
memory_min = memory_min + memory_siz
memory_max = memory_max + memory_siz
if(debug) then
write (logf,*)trim(l_here),'-d: memory_siz=',memory_siz
write (logf,*)trim(l_here),'-d: memory_min=',memory_min
write (logf,*)trim(l_here),'-d: memory_max=',memory_max
end if
! Improve
call symmetry_write
! Improve
call nuclear_potential
! Improve
call nuclear_moments
! define system here:
system%dimensions_nb = 3
call exit (l_here)
end subroutine nuclei_menu
! =========================================================================
subroutine atom_menu
! -------------------------------------------------------------------------
!
! Name : atom_menu
!
! Description : input atom
!
! Authors:X. Krokidis, F. Colonna
! Date :26 Sep 2000
! -------------------------------------------------------------------------
implicit none
! local:
character(len=max_string_len_routine_name) :: l_here
character(len=max_string_len) :: l_string
logical(bool) :: done
real(dp) :: l_real
! begin:
l_here ='atom_menu'
call enter (l_here)
! initialize:
! cur_atom = 0 in atoms_init
call check_asso_atom ('atoms', atoms, 'atoms_init', atoms_init)
! count:
cur_atom = cur_atom + 1
if(debug) write (logf,*)trim(l_here),'-d: cur_atom=',cur_atom
done = .false.
call next_command_get ('atom')
do while (.not.done)
call next_word_get ('atom')
if(trim(cur_word) .eq. 'help') then
write (logf,*) 'atom-h: atom definitions'
write (logf,*) 'atom-h: syntax:'
write (logf,*) 'atom-h: ato{m}'
write (logf,*) 'atom-h: help'
write (logf,*) 'atom-h: nam{e} = [string*8]'
write (logf,*) 'atom-h: sym{bol} = [string*2]'
write (logf,*) 'atom-h: cha{rge} = [double]'
write (logf,*) 'atom-h: mas{s} = [double]'
write (logf,*) 'atom-h: pse{udo} = [string*8]'
write (logf,*) 'atom-h: sym{bol} = [string*2] He, Al, B, ...'
write (logf,*) 'atom-h: coo{rdinates} = [list of doubles]'
write (logf,*) 'atom-h: end_ato'
write (logf,*) 'atom-h: example'
write (logf,*) 'atom-h: atom '
write (logf,*) 'atom-h: name = Nitro_1 symbol = N charge = 7 '
write (logf,*) 'atom-h: mass = 14.007 pseudo = sbk '
write (logf,*) 'atom-h: coordinates = ( 1.0, 0.5 , -3.8)'
write (logf,*) 'atom-h: end_atom'
else if(token_is('cha', l_real)) then
atoms(cur_atom)%charge = l_real
else if(token_is('coo', list_of_double)) then
call copy (atoms(cur_atom)%coordinates, list_of_double)
else if(token_is('mas', l_real)) then
atoms(cur_atom)%mass = l_real
else if(token_is('nam', l_string)) then
atoms(cur_atom)%name = string_to_upper (l_string)
else if(token_is('pse', l_string)) then
atoms(cur_atom)%pseudo_name = string_to_upper (l_string)
else if(token_is('sym', l_string)) then
atoms(cur_atom)%symbol = string_to_upper (l_string)
else
call menu_exit (l_here, done)
end if
symmetry_array(cur_atom) = cur_atom
end do
call atom_default_bynam (atoms(cur_atom)%name, cur_atom)
if(debug .or. print_level > 0) then
call atom_alone_print (cur_atom)
end if
call exit (l_here)
end subroutine atom_menu
! =========================================================================
subroutine atoms_init
! -------------------------------------------------------------------------
!
! Name : atoms_init
!
! Description : initialize atoms array
!
! Creates : atoms (atom)
!
! Authors:X. Krokidis, F. Colonna
! Date :26 Sep 2000
! -------------------------------------------------------------------------
implicit none
! local:
character(len=max_string_len_routine_name),save :: l_here
integer(i4b) :: l_atom
integer(i4b) :: atoms_nb
! header:
if (header_exe) then
l_here = "atoms_init"
! call object_create ('atoms')
return
end if
! begin:
call enter (l_here)
if(associated(atoms)) return ! improve
! count how many "atom" words are in input_words
atoms_nb = string_head_inside_cmd_count ('ato', 'nuc', 'command')
system%atoms_nb = atoms_nb
if(debug) then
write (logf,*)trim(l_here),'-d: nb of atoms found in input =',atoms_nb
end if
message = 'in routine'*l_here*' atoms_nb = 0'*CR*&
' check nuclei input'
call require (message, atoms_nb /= 0)
if (.not.associated (atoms)) then
write (logf,*)trim(l_here),': atoms not associated'
write (logf,*)trim(l_here),': allocate (atoms(',atoms_nb,')'
allocate (atoms(atoms_nb), stat = all_err)
else
write (logf,*)trim(l_here),': atoms left in memory'
endif
if (all_err /= 0) then
write (logf,'(4a)') trim(l_here), &
': allocation failed for object >atoms<'
write (logf,'(4a)') trim(l_here), &
': allocation failed for object >symmetry_array<'
write (logf,'(2a,i10)') trim(l_here), &
': dimensions are ',atoms_nb
call die (l_here, 'allocation failed')
end if
allocate (symmetry_array(atoms_nb), stat = all_err)
symmetry_type = 'default'
message = l_here+': allocate (symmetry_array('+atoms_nb+'))'
call require (message, all_err == 0)
! initialize counter:
cur_atom = 0
if(print_level > 0) write (logf,*)trim(l_here),': cur_atom intialized to 0'
! initialize atom(i)%fields
do l_atom = 1 , size(atoms)
atoms(l_atom)%pseudo_name = NOT_INIT_S
atoms(l_atom)%charge = NOT_INIT_DP
atoms(l_atom)%mass = NOT_INIT_DP
atoms(l_atom)%symbol = NOT_INIT_S
atoms(l_atom)%name = NOT_INIT_S
if(debug) write (logf,*)trim(l_here),'-d: atom # ',l_atom,' fields intialized'
end do
call exit (l_here)
end subroutine atoms_init
! =========================================================================
subroutine atom_alone_print (this_atom)
! -------------------------------------------------------------------------
!
! Name : atom_alone_print
!
! Argument (in ) : this_atom
!
! Description : prints one atom information
!
! Authors:F. Colonna
! Date :28 Feb 2001
! -------------------------------------------------------------------------
implicit none
! i/o:
integer(i4b) :: this_atom
! local:
character(len=max_string_len_routine_name) :: l_here
integer(i4b) :: i, dim_nb
real(dp) :: charge, mass , x
character(len=max_len_pseudo_name) :: pseudo
character(len=max_len_atom_name) :: name
character(len=2) :: symbol
! begin:
l_here ='atom_alone_print'
call enter (l_here)
call check_asso_atom ('atoms', atoms, 'die', die)
if(debug)write (logf,*)trim(l_here),'-d: atoms_nb =',size(atoms)
! check first:
call atom_alone_check (this_atom)
dim_nb = size(atoms(this_atom)%coordinates)
name = atoms(this_atom)%name
pseudo = atoms(this_atom)%pseudo_name
charge = atoms(this_atom)%charge
mass = atoms(this_atom)%mass
symbol = atoms(this_atom)%symbol
write (logf,'(1x,4a)') trim(l_here),': atom >',trim(name),'<'
write (logf,'(1x,4a)') trim(l_here),': pseudo >',trim(pseudo),'<'
write (logf,'(1x,4a)') trim(l_here),': symbol >',trim(symbol),'<'
write (logf,'(1x,2a,g25.15)') trim(l_here),': charge =',charge
write (logf,'(1x,2a,g25.15)') trim(l_here),': mass =',mass
system%dimensions_nb = size(atoms(1)%coordinates)
do i = 1, dim_nb
x = atoms(this_atom)%coordinates(i)
message = l_here+': coordinate # '*i*' = '*x
write (logf,'(1x,a)') trim(message)
end do
call exit (l_here)
end subroutine atom_alone_print
! =========================================================================
subroutine atoms_print
! -------------------------------------------------------------------------
!
! Name : atoms_print
!
! Needs : atoms (atom)
!
! Description : prints atoms information
!
! Authors:X. Krokidis, F. Colonna
! Date :20 Oct 2000
! Revision: 28 Feb 2001 F. Colonna
! -------------------------------------------------------------------------
implicit none
! i/o:
! local:
character(len=max_string_len_routine_name),save :: l_here
integer(i4b) :: atom_i
! header:
if (header_exe) then
l_here = "atoms_print"
call object_needed (' atoms ')
return
end if
! begin:
call enter (l_here)
call check_asso_atom ('atoms', atoms, 'die', die)
if(debug)write (logf,*)trim(l_here),'-d: atoms_nb =',size(atoms)
do atom_i = 1 , size(atoms)
call atom_alone_print (atom_i)
end do
call exit (l_here)
end subroutine atoms_print
! =========================================================================
subroutine atoms_print_res (text)
! -------------------------------------------------------------------------
!
! Name : atoms_print_res
!
! Argument (in ) : text
!
! Needs : atoms (atom)
!
! Description : prints atoms information on result file
!
! Authors:F. Colonna
! Date :26 Feb 2001
! -------------------------------------------------------------------------
implicit none
! i/o:
character(len=*) :: text
! local:
character(len=max_string_len_routine_name),save :: l_here
character(len=max_string_len) :: string
real(dp) :: charge
real(dp) :: mass
character(len=max_len_pseudo_name) :: pseudo
character(len=max_len_atom_name) :: name
character(len=2) :: symbol
integer(i4b) :: atom_i, atoms_nb
integer(i4b) :: i, dim_nb
integer(i4b) :: debut, fin
! header:
if (header_exe) then
l_here = "atoms_print_res"
call object_needed (' atoms ')
return
end if
! begin:
call enter (l_here)
include 'chksiz_atoms.inc'
atoms_nb = size (atoms)
debut = 40 - (len_trim(text)/2)
fin = 40 + (len_trim(text)/2)
if(debug) then
write (logf,*) trim(l_here),'-d: text >',trim(text),'<'
write (logf,*) trim(l_here),'-d: debut=', debut,' fin=',fin
write (logf,*) trim(l_here),'-d: atoms_nb =',size(atoms)
end if
string = dashes
string (debut-1:debut) = " "
string (debut:) = trim(text)//' '//dashes
write (res,'(a)') ' '
write (res,'(a80)') string
write (res,'(a)') ' '
! message to log file
write (logf,'(1x,2a,a80)') trim(l_here),': ',string
write(res,'(a, t15,a, t22,a, t31,a, t44,a, t58,a, t69,a, t80,a)') &
'atom# name','symbol','pseudo','charge','mass', 'x', 'y', 'z'
do atom_i = 1 , atoms_nb
name = atoms(atom_i)%name
pseudo = atoms(atom_i)%pseudo_name
charge = atoms(atom_i)%charge
mass = atoms(atom_i)%mass
symbol = atoms(atom_i)%symbol
dim_nb = size(atoms(atom_i)%coordinates)
write(res,'(i5, t7,a8, t15,a, t22,a, t26,5f11.6)') &
atom_i, name, symbol, pseudo, charge,mass, &
(atoms(atom_i)%coordinates(i),i = 1, dim_nb)
end do
call exit (l_here)
end subroutine atoms_print_res
! =========================================================================
subroutine coord_init (n_coord)
! -------------------------------------------------------------------------
!
! Name : coord_init
!
! Argument (in ) : n_coord
!
! Description : initialize coordinates array at n_coord dimension
! Remark : is not used
! Remark : this array is initialized (as it should)
! Remark : by the copy routine in atom_menu routine
!
! Authors:X. Krokidis, F. Colonna
! Date :26 Sep 2000
! -------------------------------------------------------------------------
implicit none
! i/o:
integer(i4b) :: n_coord
! local:
character(len=max_string_len_routine_name) :: l_here
! begin:
l_here ='coord_init'
call enter (l_here)
if(debug) then
write (logf,*)trim(l_here),' n_coord=',n_coord
write (logf,*)trim(l_here),' cur_atom=',cur_atom
end if
message='in routine'*l_here+' n_coord = 0'
call require (message, n_coord .ne. 0)
call allocation ('atoms%coordinates', &
atoms(cur_atom)%coordinates, n_coord)
! initialize coordinates
atoms(cur_atom)%coordinates = NOT_INIT_DP
call exit (l_here)
end subroutine coord_init
! =========================================================================
subroutine check_asso_atom (object_name, object, routine_name, routine)
! -------------------------------------------------------------------------
!
! Name : check_asso_atom
!
! Description : (object_name, object, routine_name, routine)
! Description : builds object atoms if not associated
!
! Authors:F.Colonna
! Date :21 Sep 2000
! Revision :09 Feb 2001 call die special
! -------------------------------------------------------------------------
implicit none
! i/o
character(len=*), intent (in) :: object_name, routine_name
type (atom), pointer :: object(:)
external :: routine
! local:
character(len=max_string_len_routine_name) :: l_here
! begin:
l_here ='check_asso_atom'
call enter (l_here,3)
if(debug)write (logf,*)trim(routines_21), &
'-d: object_name >',trim(object_name),'<'
if(.not.associated(object)) then
if(trim(routine_name) == "die") then
message='object >'+object_name+'< not yet created'&
+' probably because no atoms were entered'
call die (l_here, message)
else
if(debug .or. trace_level > 0) then
write (logf,'(4a)') trim(l_here), &
': object >',trim(object_name),'< not yet created'
write (logf,'(4a)') trim(l_here), &
': calling routine >',trim(routine_name),'<'
end if
call routine
end if
end if
call exit (l_here,3)
end subroutine check_asso_atom
! =========================================================================
subroutine atom_alone_check (this_atom)
! -------------------------------------------------------------------------
!
! Name : atom_alone_check
!
! Argument (in ) : this_atom
!
! Description : checks one atom information
!
! Authors:F. Colonna
! Date :28 Feb 2001
! -------------------------------------------------------------------------
implicit none
! i/o:
integer(i4b) :: this_atom
! local:
character(len=max_string_len_routine_name) :: l_here
real(dp) :: charge, mass , x
character(len=max_len_pseudo_name) :: pseudo
character(len=max_len_atom_name) :: name
character(len=2) :: symbol
integer(i4b) :: l_symb
integer(i4b) :: l_pos
integer(i4b) :: l_dim
integer(i4b) :: i
! begin:
l_here ='atom_alone_check'
call enter (l_here)
call check_asso_atom ('atoms', atoms, 'die', die)
call check_asso_mendeleiev &
('mendeleiev', mendeleiev, 'mendeleiev_init', mendeleiev_init)
if(debug)write (logf,*)trim(l_here),'-d: atoms_nb =',size(atoms)
name = atoms(this_atom)%name
pseudo = atoms(this_atom)%pseudo_name
charge = atoms(this_atom)%charge
mass = atoms(this_atom)%mass
symbol = atoms(this_atom)%symbol
if(debug) then
write (logf,*)trim(l_here),'-d: name >',trim(name),'<'
write (logf,*)trim(l_here),'-d: symbol >',trim(symbol),'<'
write (logf,*)trim(l_here),'-d: pseudo >',trim(pseudo),'<'
write (logf,*)trim(l_here),'-d: charge =',charge
write (logf,*)trim(l_here),'-d: mass =',mass
end if
! check symbol
if(symbol(1:1) == NOT_INIT_S .or. symbol(1:1) == ' ') then
message = 'in routine'*l_here*CR*'symbol not specified for atom # '&
*this_atom*' named >'+name+'<'
call require (message, symbol(1:1) /= NOT_INIT_S)
end if
! check symbol stored in table
l_symb = string_in_array_locate (symbol, mendeleiev%name, 1, n_mendeleiev)
if(debug) then
write (logf,*)trim(l_here),'-d: symb # ',l_symb
write (logf,*)trim(l_here),'-d: size(mendeleiev%name)=',size(mendeleiev%name)
do i = 1,size(mendeleiev)
write (logf,*)trim(l_here),'-d: element # ',i, &
' symbol >',trim(mendeleiev(i)%name),'<'
end do
end if
if (l_symb == 0) then
message ='in '*l_here*' atomic symbol >'*symbol*'<' &
*CR*' was not found in mendeleiev table (atom >'*name*'<)'
call require (message, l_symb /= 0)
end if
! check name
if(name(1:1) == NOT_INIT_S) then
name = symbol//'_'//encode(this_atom)
atoms(this_atom)%name = char_substitute (' ', '_', trim(name))
write (logf,*)trim(l_here),'-w: atom # ',this_atom,' name set to >',trim(name),'<'
end if
! check if name already used
if(this_atom .gt. 1) then
l_pos = string_in_array_locate (name, atoms%name, 1, this_atom-1)
if (l_pos /= 0) then
write (logf,*)trim(l_here),'-w: atom # ',this_atom, &
' has same name >',trim(name),'< as atom # ',l_pos
name = trim(name)//'_'//encode(this_atom)
atoms(this_atom)%name = char_substitute (' ', '_', trim(name))
name = atoms(this_atom)%name
write (logf,*)trim(l_here),'-w: atom # ',this_atom,' name re-set to >',trim(name),'<'
end if
end if
! check that all atoms have the same dimensions
l_dim = size(atoms(this_atom)%coordinates)
message = 'in routine'*l_here*' the number of coordinates per atom = ' &
*l_dim*' for atom # '*this_atom*CR* &
' please check coordinates input in atom menu'
call require (message, l_dim > 0)
message = 'in routine'*l_here*' the number of coordinates per atom = ' &
*l_dim*' for atom # '*this_atom*' and ='* &
size(atoms(1)%coordinates)*' for atom # 1'
call require (message, l_dim == size(atoms(1)%coordinates) )
! check x
do i = 1,size(atoms(this_atom)%coordinates)
x = atoms(this_atom)%coordinates(i)
! require_x:
if(x <= NOT_INIT_DP) then
message = ' atom #'*this_atom*' >'*name* &
'< has uninitialized coordinate # '*i*'='*trim(double_to_string(x))
call die (l_here, message)
end if
enddo
! check mass
if(mass <= NOT_INIT_DP) then
message = 'mass not specified for atom # '*this_atom* &
' >'+name+'< taken from mendeleiev table'
atoms(this_atom)%mass = mendeleiev(l_symb)%mass
if(warning_level .ge. 1) call warning (l_here, message)
end if
! check charge
if(charge <= NOT_INIT_DP) then
message = 'charge not specified for atom # '*this_atom* &
' >'+name+'< taken from mendeleiev table'
atoms(this_atom)%charge = mendeleiev(l_symb)%zed
if(warning_level .ge. 1) call warning (l_here, message)
end if
call exit (l_here)
end subroutine atom_alone_check
! =========================================================================
subroutine atoms_check
! -------------------------------------------------------------------------
!
! Name : atoms_check
!
! Description : checks that all atoms(i)%fields are set
!
! Authors:X. Krokidis, F. Colonna
! Date :20 Oct 2000
! -------------------------------------------------------------------------
implicit none
! i/o:
! local:
character(len=max_string_len_routine_name) :: l_here
integer(i4b) :: l_atom
! begin:
l_here ='atoms_check'
call enter (l_here)
call check_asso_mendeleiev &
('mendeleiev', mendeleiev, 'mendeleiev_init', mendeleiev_init)
call check_asso_atom ('atoms', atoms, 'die', die)
do l_atom = 1 , size(atoms)
call atom_alone_check (l_atom)
end do
call exit (l_here)
end subroutine atoms_check
! =========================================================================
subroutine nuclear_potential
! -------------------------------------------------------------------------
!
! Name : nuclear_potential
!
! Description : calculates the nuclear potential
!
! Authors:X. Krokidis, P. Pernot
! Date :23 Jan 2001
! Revision: 26 Feb 2001 F.Colonna call check_asso_atom ....
! -------------------------------------------------------------------------
implicit none
! i/o:
! local:
character(len=max_string_len_routine_name) :: l_here
integer(i4b) :: i, j, idim
integer(i4b) :: dimensions_nb
real(dp) :: potential , sum_of_mass
real(dp), pointer :: R(:,:) => NULL()
! begin:
l_here ='nuclear_potential'
call enter (l_here)
call check_asso_atom ('atoms', atoms, 'atoms_init', atoms_init)
if(debug) write (logf,*)trim(l_here),' atoms_nb=',size(atoms)
dimensions_nb = size(atoms(1)%coordinates)
if(system%dimensions_nb > 0 .and. system%dimensions_nb /= dimensions_nb) then
call die (l_here, ' incorrect number of dimensions')
else
system%dimensions_nb = size(atoms(1)%coordinates)
end if
if(debug) write (logf,*)trim(l_here),' dimensions_nb=',dimensions_nb
! computation of the center of mass
message = 'in routine'*l_here*' dimensions_nb = '+dimensions_nb
call require (message, dimensions_nb > 0)
allocate (center_of_mass(dimensions_nb), stat = all_err)
message = l_here+': allocate (center_of_mass('+dimensions_nb+'))'
call require (message, all_err == 0)
if(print_level > 1) write (logf,*)trim(message)
center_of_mass(:) = ZERO
sum_of_mass = ZERO
do i = 1, size(atoms)
do j = 1,dimensions_nb
center_of_mass(j) = center_of_mass(j) + atoms(i)%mass * atoms(i)%coordinates(j)
enddo
sum_of_mass = sum_of_mass + atoms(i)%mass
enddo
center_of_mass = center_of_mass/sum_of_mass
! write (logf,'(a,1g15.8)')'center of mass. X = ',center_of_mass(1)
! write (logf,'(a,1g15.8)')' Y = ',center_of_mass(2)
! write (logf,'(a,1g15.8)')' Z = ',center_of_mass(3)
allocate(R(size(atoms),size(atoms)))
potential = ZERO
do i = 1, size(atoms)-1
do j = i+1, size(atoms)
R(i, j) = ZERO
do idim = 1, 3
R(i, j) = R(i, j) + (atoms(i)%coordinates(idim)-atoms(j)%coordinates(idim))**2
enddo
R(i,j) = sqrt(R(i,j))
potential = potential + atoms(i)%charge * atoms(j)%charge / R(i,j)
enddo
enddo
if (size(atoms) /= 1) then
! write (logf,'(a)')'Nuclear distances matrix :'
! write (logf,'(5X,10(3x,a2,3x))') (atoms(i)%symbol,i=1,size(atoms))
! do i=1,size(atoms)
! write (logf,'(a2,3x,10f8.4)') atoms(i)%symbol,(R(j,i),j=1,i-1),0.
! enddo
! write (logf,'(a,g15.8)')'Nuclear potential energy = ',potential
else
potential = ZERO
endif
nuclear_potential_energy = potential
call exit (l_here)
end subroutine nuclear_potential
! =========================================================================
subroutine nuclear_moments
! -------------------------------------------------------------------------
!
! Name : nuclear_moments
!
! Description : calculates the nuclear dipole moment
! Description : calculates the nuclear quadrupole moment
!
! Authors:X. Krokidis
! Date :07 Feb 2001
! -------------------------------------------------------------------------
implicit none
! i/o:
! local:
character(len=max_string_len_routine_name) :: l_here
integer(i4b) :: dimensions_nb
integer(i4b) :: i, j
! begin:
l_here ='nuclear_moments'
call enter (l_here)
call check_asso_atom ('atoms', atoms, 'atoms_init', atoms_init)
dimensions_nb = size(atoms(1)%coordinates)
allocate (nuclear_dipole_moment(dimensions_nb))
nuclear_dipole_moment = ZERO
do j = 1,dimensions_nb
do i = 1, size(atoms)
nuclear_dipole_moment(j) = nuclear_dipole_moment(j) &
+ atoms(i)%charge &
*(atoms(i)%coordinates(j)-center_of_mass(j))
enddo
enddo
write (logf,'(1x,2a,3(a,g13.6))') trim(l_here), &
': Nuclear dipole moment :', &
' x=',nuclear_dipole_moment(1), &
' y=',nuclear_dipole_moment(2), &
' z=',nuclear_dipole_moment(3)
call exit (l_here)
end subroutine nuclear_moments
! =========================================================================
subroutine symmetry_write
! -------------------------------------------------------------------------
!
! Name : symmetry_write
!
! Description : writes symmetry in output file
!
! Authors:X. Krokidis
! Date :23 Jan 2001
! -------------------------------------------------------------------------
implicit none
! i/o:
! local:
character(len=max_string_len_routine_name) :: l_here
! begin:
l_here ='symmetry_write'
call enter (l_here)
if(debug) write (logf,*)trim(l_here),'-d: '//trim(symmetry_type)//' symmetry entered'
call printarray (l_here, symmetry_array, 'symmetry correspondence:')
call exit (l_here)
end subroutine symmetry_write
! =========================================================================
subroutine atom_default_bynam (name, atom_n)
! -------------------------------------------------------------------------
!
! Name : atom_default_bynam
!
! Argument (in ) : name
! Argument (in ) : atom_n atom number from menu
!
! Description : fill atom array element with default from mendeleiev
!
! Authors:F. Colonna
! Date :28 Jan 2003
! -------------------------------------------------------------------------
implicit none
! i/o:
character(len=*) :: name
! local:
character(len=max_string_len_routine_name),save :: l_here
character(len=max_len_atom_symbol) :: symb
character(len=max_len_atom_name) :: temp
integer(i4b) :: atom_i, atoms_nb, atom_n
integer(i4b) :: symb_i, symbs_nb
! begin:
l_here = "atom_default_bynam"
call enter (l_here)
if(debug) then
write (logf,*)trim(l_here),'-d: atom name >',trim(name),'<'
write (logf,*)trim(l_here),'-d: atom number =',atom_n
end if
atoms_nb = size (atoms)
atom_i = string_in_array_locate_or_die (name , atoms%name, 1, atoms_nb)
if (atom_i /= atom_n) then
temp = name
name = temp+"_"+atom_n
message = 'atom # '*atom_i*' and atom # '*atom_n &
*' have same name >'+temp+'< changed to >'+name+'<'
call warning (l_here, message)
end if
call check_asso_mendeleiev &
('mendeleiev', mendeleiev, 'mendeleiev_init', mendeleiev_init)
symbs_nb = size (mendeleiev)
if (atoms(atom_n)%symbol == NOT_INIT_S ) then
symb = name(1:max_len_atom_symbol)
else
symb = atoms(atom_n)%symbol
end if
symb_i = string_in_array_locate (symb , mendeleiev%name, 1, symbs_nb)
if (symb_i == 0) then
message = 'in routine'*l_here*' symbol '*symb*CR*&
'not in mendeleiev table'
call require (message, symb_i /= 0)
else
atoms(atom_n)%symbol = symb
end if
if(debug) then
write (logf,*)trim(l_here),'-d: atom # ',atom_n
write (logf,*)trim(l_here),'-d: atom name >',trim(name),'<'
write (logf,*)trim(l_here),'-d: symbol # ',symb_i
write (logf,*)trim(l_here),'-d: atom symbol >',trim(symb),'<'
end if
if (atoms(atom_n)%mass <= NOT_INIT_DP ) then
atoms(atom_n)%mass = mendeleiev(symb_i)%mass
end if
if (atoms(atom_n)%charge <= NOT_INIT_DP ) then
atoms(atom_n)%charge = mendeleiev(symb_i)%zed
end if
call exit (l_here)
end subroutine atom_default_bynam
! ==========================================================================
function atom_close_to_point (x, y, z) result (result)
! --------------------------------------------------------------------------
!
! Name : atom_close_to_point
!
! Argument (in ) : x, y, z
!
! Descriptions : index of an atom at less than tiny(1.d0) from (x,y,z)
! Descriptions : zero is returned if no index found
!
! Authors : F. Colonna
! Date : 14 Mar 2003
! --------------------------------------------------------------------------
implicit none
! result:
integer(i4b) :: result
! i/o:
real(dp) :: x, y, z
! local:
character(len=max_string_len_routine_name),save :: l_here
integer(i4b) :: atom_i, atoms_nb
integer(i4b) :: dimensions_nb
real(dp) :: dx, dy, dz
real(dp) :: dist2
! begin:
l_here = "atom_close_to_point"
call enter (l_here)
atoms_nb = size (atoms)
dimensions_nb = size (atoms(1)%coordinates)
if(atoms_nb == 0) then
if(atoms_nb <= 0) call die (l_here, ' array of coordinates ='+atoms_nb)
end if
if (dimensions_nb /= 3) then
call die (l_here, ' 3D space is required for this type of system')
end if
result = 0
do atom_i = 1, atoms_nb
dx = x - atoms(atom_i)%coordinates (1)
dy = y - atoms(atom_i)%coordinates (2)
dz = z - atoms(atom_i)%coordinates (3)
dist2 = dx*dx + dy*dy + dz*dz
if (dist2 .lt. tiny(1.d0)) then
result = atom_i
exit
end if
end do
if(debug) then
write (logf,*)trim(l_here),'-d: result =',result
end if
call exit (l_here)
end function atom_close_to_point
end module nuclei_mod
back to top