From e887d2c5ab82154dc0146e0cbeaa9f188aefe88b Mon Sep 17 00:00:00 2001 From: Dynamics of Condensed Matter <30792324+DCM-Uni-Paderborn@users.noreply.github.com> Date: Sat, 25 Apr 2026 15:11:05 +0200 Subject: [PATCH] FIST: add silicon EIP Stillinger-Weber and Tersoff models (#5095) Co-authored-by: Thomas D. Kuehne --- src/common/bibliography.F | 34 +- src/eip_environment.F | 8 +- src/eip_environment_types.F | 14 +- src/eip_silicon.F | 2340 ++++++++++++++++- src/force_env_methods.F | 20 +- src/input_constants.F | 4 +- src/input_cp2k_eip.F | 37 +- tests/Fist/regtest-1-1/Si_1000_bazant.inp | 30 + tests/Fist/regtest-1-1/Si_1000_lenosky.inp | 30 + .../regtest-1-1/Si_1000_stillinger_weber.inp | 30 + tests/Fist/regtest-1-1/Si_1000_tersoff.inp | 30 + tests/Fist/regtest-1-1/TEST_FILES.toml | 5 +- .../Si_1000.inp => sample_xyz/Si_1000.xyz} | 30 +- 13 files changed, 2544 insertions(+), 68 deletions(-) create mode 100644 tests/Fist/regtest-1-1/Si_1000_bazant.inp create mode 100644 tests/Fist/regtest-1-1/Si_1000_lenosky.inp create mode 100644 tests/Fist/regtest-1-1/Si_1000_stillinger_weber.inp create mode 100644 tests/Fist/regtest-1-1/Si_1000_tersoff.inp rename tests/Fist/{regtest-1-1/Si_1000.inp => sample_xyz/Si_1000.xyz} (99%) diff --git a/src/common/bibliography.F b/src/common/bibliography.F index f2707eaaa6..e5957772f5 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -44,7 +44,8 @@ MODULE bibliography VandeVondele2007, Ortiz1994, Becke1988, Perdew1996, Zhang1998, & Perdew2008, Lee1988, Heyd2006, Vydrov2006, Heyd2003, Heyd2004, & Vosko1980, Aguado2003, Essmann1995, Ewald1921, Darden1993, & - Siepmann1995, Tersoff1988, Tosi1964a, Tosi1964b, Yamada2000, & + Bazant1996, Bazant1997, Goedecker2002, Lenosky2000, Siepmann1995, & + Stillinger1985, Tersoff1988, Tosi1964a, Tosi1964b, Yamada2000, & Dudarev1997, Dudarev1998, Dewar1977, Dewar1985, Rocha2006, & Stewart1989, Thiel1992, Repasky2002, Stewart2007, Weber2008, & Hunt2003, Guidon2008, Elber1987, Jonsson1998, Jonsson2000_1, & @@ -327,12 +328,43 @@ CONTAINS source="J. Phys. Chem. Solids", volume="25", pages="45-52", & year=1964, doi="10.1016/0022-3697(64)90160-X") + CALL add_reference(key=Stillinger1985, & + authors=s2a("F. H. Stillinger", "T. A. Weber"), & + title="Computer simulation of local order in condensed phases of silicon", & + source="Phys. Rev. B", volume="31", pages="5262-5271", & + year=1985, doi="10.1103/PhysRevB.31.5262") + CALL add_reference(key=Tersoff1988, & authors=s2a("J. Tersoff"), & title="Empirical interatomic potential for silicon with improved elastic properties", & source="Phys. Rev. B", volume="38", pages="9902-9905", & year=1988, doi="10.1103/PhysRevB.38.9902") + CALL add_reference(key=Bazant1996, & + authors=s2a("M. Z. Bazant", "E. Kaxiras"), & + title="Modeling of covalent bonding in solids by inversion of cohesive energy curves", & + source="Phys. Rev. Lett.", volume="77", pages="4370-4373", & + year=1996, doi="10.1103/PhysRevLett.77.4370") + + CALL add_reference(key=Bazant1997, & + authors=s2a("M. Z. Bazant", "E. Kaxiras", "J. F. Justo"), & + title="Environment-dependent interatomic potential for bulk silicon", & + source="Phys. Rev. B", volume="56", pages="8542-8552", & + year=1997, doi="10.1103/PhysRevB.56.8542") + + CALL add_reference(key=Lenosky2000, & + authors=s2a("T. J. Lenosky", "B. Sadigh", "E. Alonso", "V. V. Bulatov", & + "T. Diaz de la Rubia", "J. Kim", "A. F. Voter", "J. D. Kress"), & + title="Highly optimized empirical potential model of silicon", & + source="Model. Simul. Mater. Sci. Eng.", volume="8", pages="825-841", & + year=2000, doi="10.1088/0965-0393/8/6/305") + + CALL add_reference(key=Goedecker2002, & + authors=s2a("S. Goedecker"), & + title="Optimization and parallelization of a force field for silicon using OpenMP", & + source="Comput. Phys. Commun.", volume="148", pages="124-135", & + year=2002, doi="10.1016/S0010-4655(02)00466-6") + CALL add_reference(key=Siepmann1995, & authors=s2a("J. I. Siepmann", "M. Sprik"), & title="Influence of surface topology and electrostatic potential on water/electrode systems ", & diff --git a/src/eip_environment.F b/src/eip_environment.F index 5cbb7ae8b5..a09226d650 100644 --- a/src/eip_environment.F +++ b/src/eip_environment.F @@ -9,7 +9,7 @@ !> \brief Methods and functions on the EIP environment !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** MODULE eip_environment USE atomic_kind_types, ONLY: atomic_kind_type,& @@ -60,7 +60,7 @@ CONTAINS !> \param subsys_section ... !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_init(eip_env, root_section, para_env, force_env_section, & subsys_section) @@ -107,7 +107,7 @@ CONTAINS !> \param subsys_section ... !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_init_subsys(eip_env, subsys, subsys_section) TYPE(eip_environment_type), POINTER :: eip_env @@ -186,7 +186,7 @@ CONTAINS !> \param eip_env The eip environment to retain !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_init_model(eip_env) TYPE(eip_environment_type), POINTER :: eip_env diff --git a/src/eip_environment_types.F b/src/eip_environment_types.F index 24de373702..72d4cc73be 100644 --- a/src/eip_environment_types.F +++ b/src/eip_environment_types.F @@ -9,7 +9,7 @@ !> \brief The environment for the empirical interatomic potential methods. !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** MODULE eip_environment_types USE atomic_kind_list_types, ONLY: atomic_kind_list_create,& @@ -78,7 +78,7 @@ MODULE eip_environment_types !> \param virial Dummy virial pointer !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** TYPE eip_environment_type INTEGER :: eip_model = 0 @@ -104,7 +104,7 @@ CONTAINS !> \param eip_env The eip environment to release !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_env_release(eip_env) @@ -162,7 +162,7 @@ CONTAINS !> eip_environment_type !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_env_get(eip_env, eip_model, eip_energy, eip_energy_var, & eip_forces, coord_avg, coord_var, count, subsys, & @@ -265,7 +265,7 @@ CONTAINS !> \param eip_potential_energy The EIP potential energy !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) !> \note !> For possible missing arguments see the attributes of eip_environment_type ! ************************************************************************************************** @@ -367,7 +367,7 @@ CONTAINS !> \param eip_env The eip environment to be reinitialized !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_env_clear(eip_env) @@ -403,7 +403,7 @@ CONTAINS !> \param eip_env The eip environment to be created !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_env_create(eip_env) diff --git a/src/eip_silicon.F b/src/eip_silicon.F index e723c40e9a..1e5bf4d308 100644 --- a/src/eip_silicon.F +++ b/src/eip_silicon.F @@ -12,7 +12,7 @@ !> empirical interatomic potentials for Silicon. !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** MODULE eip_silicon USE atomic_kind_list_types, ONLY: atomic_kind_list_type @@ -48,7 +48,7 @@ MODULE eip_silicon CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eip_silicon' ! *** Public subroutines *** - PUBLIC :: eip_bazant, eip_lenosky + PUBLIC :: eip_bazant, eip_lenosky, eip_stillinger_weber, eip_tersoff !*** @@ -69,7 +69,7 @@ CONTAINS !> using OpenMP; CPC 148, 1 (2002) !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_bazant(eip_env) TYPE(eip_environment_type), POINTER :: eip_env @@ -235,7 +235,7 @@ CONTAINS !> using OpenMP; CPC 148, 1 (2002) !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_lenosky(eip_env) TYPE(eip_environment_type), POINTER :: eip_env @@ -391,13 +391,334 @@ CONTAINS END SUBROUTINE eip_lenosky +! ************************************************************************************************** +!> \brief Interface routine of the Stillinger-Weber force field to CP2K +!> \param eip_env ... +!> \par Literature +!> F.H. Stillinger and T.A. Weber: +!> Computer simulation of local order in condensed phases of silicon; +!> Phys. Rev. B 31, 5262 (1985) +!> \par History +!> 04.2026 added [Thomas D. Kuehne, tkuehne@cp2k.org] +! ************************************************************************************************** + SUBROUTINE eip_stillinger_weber(eip_env) + TYPE(eip_environment_type), POINTER :: eip_env + + CHARACTER(len=*), PARAMETER :: routineN = 'eip_stillinger_weber' + + INTEGER :: handle, i, iparticle, iparticle_kind, & + iparticle_local, iw, natom, & + nparticle_kind, nparticle_local + REAL(KIND=dp) :: ekin, ener, mass + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz + REAL(KIND=dp), DIMENSION(3) :: abc + TYPE(atomic_kind_list_type), POINTER :: atomic_kinds + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(distribution_1d_type), POINTER :: local_particles + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(section_vals_type), POINTER :: eip_section + + CALL timeset(routineN, handle) + + NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, & + atomic_kind, local_particles, subsys, atomic_kind_set, para_env) + + ekin = 0.0_dp + + logger => cp_get_default_logger() + + CPASSERT(ASSOCIATED(eip_env)) + + CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, & + subsys=subsys, local_particles=local_particles, & + atomic_kind_set=atomic_kind_set) + CALL get_cell(cell=cell, abc=abc) + + eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP") + natom = SIZE(particle_set) + + ALLOCATE (rxyz(3, natom)) + + DO i = 1, natom + rxyz(:, i) = particle_set(i)%r(:)*angstrom + END DO + + CALL eip_stillinger_weber_silicon(nat=natom, alat=abc*angstrom, & + rxyz0=rxyz, fxyz=eip_env%eip_forces, & + etot=ener, count=eip_env%count) + + eip_env%coord_avg = 0.0_dp + eip_env%coord_var = 0.0_dp + + CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds) + + nparticle_kind = atomic_kinds%n_els + + DO iparticle_kind = 1, nparticle_kind + atomic_kind => atomic_kind_set(iparticle_kind) + CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) + nparticle_local = local_particles%n_el(iparticle_kind) + DO iparticle_local = 1, nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + ekin = ekin + 0.5_dp*mass* & + (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) & + + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) & + + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3)) + END DO + END DO + + CALL cp_subsys_get(subsys=subsys, para_env=para_env) + CALL para_env%sum(ekin) + eip_env%eip_kinetic_energy = ekin + + eip_env%eip_potential_energy = ener/evolt + eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy + eip_env%eip_energy_var = 0.0_dp + + DO i = 1, natom + particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom + END DO + + DEALLOCATE (rxyz) + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", & + extension=".mmLog") + + CALL eip_print_energies(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%ENERGIES") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", & + extension=".mmLog") + + CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%ENERGIES_VAR") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%FORCES"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", & + extension=".mmLog") + + CALL eip_print_forces(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%FORCES") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", & + extension=".mmLog") + + CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COORD_AVG") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", & + extension=".mmLog") + + CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COORD_VAR") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COUNT"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", & + extension=".mmLog") + + CALL eip_print_count(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COUNT") + END IF + + CALL timestop(handle) + + END SUBROUTINE eip_stillinger_weber + +! ************************************************************************************************** +!> \brief Interface routine of the Tersoff force field to CP2K +!> \param eip_env ... +!> \par Literature +!> J. Tersoff: +!> New empirical approach for the structure and energy of covalent systems; +!> Phys. Rev. Lett. 61, 2879 (1988) +!> J. Tersoff: +!> Modeling solid-state chemistry: Interatomic potentials for multicomponent systems; +!> Phys. Rev. B 39, 5566 (1989) +!> \par History +!> 04.2026 added [Thomas D. Kuehne, tkuehne@cp2k.org] +! ************************************************************************************************** + SUBROUTINE eip_tersoff(eip_env) + TYPE(eip_environment_type), POINTER :: eip_env + + CHARACTER(len=*), PARAMETER :: routineN = 'eip_tersoff' + + INTEGER :: handle, i, iparticle, iparticle_kind, & + iparticle_local, iw, natom, & + nparticle_kind, nparticle_local + REAL(KIND=dp) :: ekin, ener, mass + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz + REAL(KIND=dp), DIMENSION(3) :: abc + TYPE(atomic_kind_list_type), POINTER :: atomic_kinds + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(distribution_1d_type), POINTER :: local_particles + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(section_vals_type), POINTER :: eip_section + + CALL timeset(routineN, handle) + + NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, & + atomic_kind, local_particles, subsys, atomic_kind_set, para_env) + + ekin = 0.0_dp + + logger => cp_get_default_logger() + + CPASSERT(ASSOCIATED(eip_env)) + + CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, & + subsys=subsys, local_particles=local_particles, & + atomic_kind_set=atomic_kind_set) + CALL get_cell(cell=cell, abc=abc) + + eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP") + natom = SIZE(particle_set) + + ALLOCATE (rxyz(3, natom)) + + DO i = 1, natom + rxyz(:, i) = particle_set(i)%r(:)*angstrom + END DO + + CALL eip_tersoff_silicon(nat=natom, alat=abc*angstrom, rxyz=rxyz, & + fxyz=eip_env%eip_forces, etot=ener, & + count=eip_env%count) + + eip_env%coord_avg = 0.0_dp + eip_env%coord_var = 0.0_dp + + CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds) + + nparticle_kind = atomic_kinds%n_els + + DO iparticle_kind = 1, nparticle_kind + atomic_kind => atomic_kind_set(iparticle_kind) + CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) + nparticle_local = local_particles%n_el(iparticle_kind) + DO iparticle_local = 1, nparticle_local + iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) + ekin = ekin + 0.5_dp*mass* & + (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) & + + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) & + + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3)) + END DO + END DO + + CALL cp_subsys_get(subsys=subsys, para_env=para_env) + CALL para_env%sum(ekin) + eip_env%eip_kinetic_energy = ekin + + eip_env%eip_potential_energy = ener/evolt + eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy + eip_env%eip_energy_var = 0.0_dp + + DO i = 1, natom + particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom + END DO + + DEALLOCATE (rxyz) + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", & + extension=".mmLog") + + CALL eip_print_energies(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%ENERGIES") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", & + extension=".mmLog") + + CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%ENERGIES_VAR") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%FORCES"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", & + extension=".mmLog") + + CALL eip_print_forces(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%FORCES") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", & + extension=".mmLog") + + CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COORD_AVG") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", & + extension=".mmLog") + + CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COORD_VAR") + END IF + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + eip_section, "PRINT%COUNT"), cp_p_file)) THEN + iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", & + extension=".mmLog") + + CALL eip_print_count(eip_env=eip_env, output_unit=iw) + CALL cp_print_key_finished_output(iw, logger, eip_section, & + "PRINT%COUNT") + END IF + + CALL timestop(handle) + + END SUBROUTINE eip_tersoff + ! ************************************************************************************************** !> \brief Print routine for the EIP energies !> \param eip_env The eip environment of matter !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) !> \note !> As usual the EIP energies differ from the DFT energies! !> Only the relative energy differences are correctly reproduced. @@ -423,7 +744,7 @@ CONTAINS !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_print_energy_var(eip_env, output_unit) TYPE(eip_environment_type), POINTER :: eip_env @@ -452,7 +773,7 @@ CONTAINS !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_print_forces(eip_env, output_unit) TYPE(eip_environment_type), POINTER :: eip_env @@ -491,7 +812,7 @@ CONTAINS !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_print_coord_avg(eip_env, output_unit) TYPE(eip_environment_type), POINTER :: eip_env @@ -520,7 +841,7 @@ CONTAINS !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_print_coord_var(eip_env, output_unit) TYPE(eip_environment_type), POINTER :: eip_env @@ -549,7 +870,7 @@ CONTAINS !> \param output_unit The output unit !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_print_count(eip_env, output_unit) TYPE(eip_environment_type), POINTER :: eip_env @@ -599,7 +920,7 @@ CONTAINS !> using OpenMP; CPC 148, 1 (2002) !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_bazant_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, & coord_var, count) @@ -1764,7 +2085,7 @@ CONTAINS !> using OpenMP; CPC 148, 1 (2002) !> \par History !> 03.2006 initial create [tdk] -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE eip_lenosky_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, & coord_var, count) @@ -3000,4 +3321,1999 @@ CONTAINS RETURN END SUBROUTINE splint +! ************************************************************************************************** +! Additional EIP kernels consolidated here to match the original Bazant/Lenosky layout. +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param nat ... +!> \param alat ... +!> \param rxyz0 ... +!> \param fxyz ... +!> \param etot ... +!> \param count ... +! ************************************************************************************************** + SUBROUTINE eip_stillinger_weber_silicon(nat, alat, rxyz0, fxyz, etot, count) +!***************************************************************************************** +! This subroutine evaluates the Stillinger Weber Silicon potential with linear scaling +! COPYRIGHT +! Copyright (C) 2009 AIST, UNIBAS +! This file is distributed under the terms of the +! GNU General Public License, see +! http://www.gnu.org/copyleft/gpl.txt . +! +! Implementation: Original version was written by Tetsuya Morishita, AIST Tsukuba (JP) +! Improved by M. Amsler, S. Goedecker, Basel University (CH), 2009 +! +! Note: +! +! aa is the parameter A given on page 5263 of PRB 31, 5262 (1985). +! bb is the parameter B given on page 5263 of PRB 31, 5262 (1985). +! ra is the parameter a given on page 5263 of PRB 31, 5262 (1985). +! gam and ramda are the parameters gamma and lambda +! for the 3-body term, respectively (see Eq. (2.5) in the paper). +! +! Input: +! nat, integer: the number of atoms +! alat, real(8), dim(3) : the three edges of the orthoromic simulation cell, periodic boundaries are applied +! and atoms outside the cell will be brought back into the box +! rxyz, real(8), dim(3,nat) : the xyz cartesian components of the atomic positions in Angstroem +! +! Output: +! fxyz, real(8), dim(3,nat): the xyz cartesian forces in on the corresponding atomic components n eV/A +! etot, real(8) : total potential energy, 2-body and 3-body, in eV +! count, real(8) : increased by 1.d0 at each call of this subroutine, +! needs to be initialized to 0.d0 before calling this routine for the first time +! +! Other variables: +! p: the 2-body potential energy +! p3: the 3-body potential energy +! fx(i) , fy(i) , fz(i) are the 2-body forces on atom i. +! fx3(i), fy3(i), fz3(i) are the 3-body forces on atom i. +! fxyz(3,nat) contains both 2-body and 3-body forces +! +! All units follow the description in PRB 31, 5262 (1985). +!***************************************************************************************** + + INTEGER :: nat + REAL(8) :: alat(3), rxyz0(3, nat), fxyz(3, nat), & + etot, count + + INTEGER :: i, iam, iat, ii, il, in, indlst, & + indlstx, ipb, l1, l2, l3, laymx, ll1, & + ll2, ll3, myspace, myspaceout, ncx, & + ndat, nn, nnbrx, npjkx, npjx, npr + INTEGER, ALLOCATABLE, DIMENSION(:) :: lay, lstb + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta + INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell + REAL(8) :: cut, cut2, eps, esigma, fx(nat), & + fx3(nat), fy(nat), fy3(nat), fz(nat), & + fz3(nat), isigma, p, p3, pv3, ra, & + rlc1i, rlc2i, rlc3i, sigma + REAL(8), ALLOCATABLE, DIMENSION(:, :) :: rel, rxyz + + PARAMETER(ra=1.8d0) + PARAMETER(sigma=2.0951d0, eps=2.167239428587d0) + + count = count + 1.d0 + cut = sigma*ra*2.d0 + isigma = 1.d0/sigma + esigma = eps*isigma + +! linear scaling calculation of verlet list, only serial + ll1 = INT(alat(1)/cut) + IF (ll1 < 1) CPABORT("alat(1) too small") + ll2 = INT(alat(2)/cut) + IF (ll2 < 1) CPABORT("alat(2) too small") + ll3 = INT(alat(3)/cut) + IF (ll3 < 1) CPABORT("alat(3) too small") + + npr = 1 + ncx = 29 + DO + ncx = ncx*2 + ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3)) + icell(0, :, :, :) = 0 + rlc1i = ll1/alat(1) + rlc2i = ll2/alat(2) + rlc3i = ll3/alat(3) + + DO iat = 1, nat + rxyz0(1, iat) = MODULO(MODULO(rxyz0(1, iat), alat(1)), alat(1)) + rxyz0(2, iat) = MODULO(MODULO(rxyz0(2, iat), alat(2)), alat(2)) + rxyz0(3, iat) = MODULO(MODULO(rxyz0(3, iat), alat(3)), alat(3)) + l1 = INT(rxyz0(1, iat)*rlc1i) + l2 = INT(rxyz0(2, iat)*rlc2i) + l3 = INT(rxyz0(3, iat)*rlc3i) + + ii = icell(0, l1, l2, l3) + ii = ii + 1 + icell(0, l1, l2, l3) = ii + IF (ii > ncx) THEN + DEALLOCATE (icell) + EXIT + END IF + icell(ii, l1, l2, l3) = iat + END DO + IF (ALLOCATED(icell)) EXIT + END DO + +! duplicate all atoms within boundary layer + laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8) + nn = nat + laymx + ALLOCATE (rxyz(3, nn), lay(nn)) + DO iat = 1, nat + lay(iat) = iat + rxyz(1, iat) = rxyz0(1, iat) + rxyz(2, iat) = rxyz0(2, iat) + rxyz(3, iat) = rxyz0(3, iat) + END DO + il = nat +! xy plane + DO l2 = 0, ll2 - 1 + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, l2, 0) + icell(0, l1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, l1, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, l2, ll3 - 1) + icell(0, l1, l2, -1) = in + DO ii = 1, in + i = icell(ii, l1, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + END DO + +! yz plane + DO l3 = 0, ll3 - 1 + DO l2 = 0, ll2 - 1 + + in = icell(0, 0, l2, l3) + icell(0, ll1, l2, l3) = in + DO ii = 1, in + i = icell(ii, 0, l2, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, l2, l3) + icell(0, -1, l2, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + END DO + +! xz plane + DO l3 = 0, ll3 - 1 + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, 0, l3) + icell(0, l1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, l1, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, l1, ll2 - 1, l3) + icell(0, l1, -1, l3) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + END DO + +! x axis + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, 0, 0) + icell(0, l1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, l1, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, 0, ll3 - 1) + icell(0, l1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, l1, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, l1, ll2 - 1, 0) + icell(0, l1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, ll2 - 1, ll3 - 1) + icell(0, l1, -1, -1) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + +! y axis + DO l2 = 0, ll2 - 1 + + in = icell(0, 0, l2, 0) + icell(0, ll1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, 0, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, l2, ll3 - 1) + icell(0, ll1, l2, -1) = in + DO ii = 1, in + i = icell(ii, 0, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, l2, 0) + icell(0, -1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, l2, ll3 - 1) + icell(0, -1, l2, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + +! z axis + DO l3 = 0, ll3 - 1 + + in = icell(0, 0, 0, l3) + icell(0, ll1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, 0, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, 0, l3) + icell(0, -1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, 0, ll2 - 1, l3) + icell(0, ll1, -1, l3) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, l3) + icell(0, -1, -1, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + +! corners + in = icell(0, 0, 0, 0) + icell(0, ll1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, 0, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, 0, 0) + icell(0, -1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, ll2 - 1, 0) + icell(0, ll1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, 0) + icell(0, -1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, 0, ll3 - 1) + icell(0, ll1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, 0, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, 0, ll3 - 1) + icell(0, -1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, 0, ll2 - 1, ll3 - 1) + icell(0, ll1, -1, -1) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1) + icell(0, -1, -1, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + ALLOCATE (lsta(2, nat)) + nnbrx = 300 + DO + nnbrx = 3*nnbrx/2 + ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat)) + + indlstx = 0 + + npr = 1 + iam = 0 + + cut2 = cut**2 +! assign contiguous portions of the arrays lstb and rel to the threads (this version only contains one thread) + myspace = (nat*nnbrx)/npr + IF (iam == 0) myspaceout = myspace +! Verlet list, relative positions + indlst = 0 + DO l3 = 0, ll3 - 1 + DO l2 = 0, ll2 - 1 + DO l1 = 0, ll1 - 1 + DO ii = 1, icell(0, l1, l2, l3) + iat = icell(ii, l1, l2, l3) + IF (((iat - 1)*npr)/nat == iam) THEN + lsta(1, iat) = iam*myspace + indlst + 1 + CALL sw_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, & + rxyz, icell, lstb(iam*myspace + 1), lay, rel(1, iam*myspace + 1), cut2, indlst) + lsta(2, iat) = iam*myspace + indlst + ipb = lsta(1, iat) + ndat = lsta(2, iat) - lsta(1, iat) + 1 + END IF + END DO + END DO + END DO + END DO + indlstx = MAX(indlstx, indlst) + + IF (indlstx < myspaceout) EXIT + DEALLOCATE (lstb, rel) + END DO + + npr = 1 + iam = 0 + npjx = 300; npjkx = 6000 +!end of creating pairlist part------------------------------------------------------------ + +!start energy and force calculation------------------------------------------------------- +!set all variables to zero + p = 0.0d0 + p3 = 0.0d0 + pv3 = 0.0d0 + fx(:) = 0.0d0 + fy(:) = 0.0d0 + fz(:) = 0.0d0 + fx3(:) = 0.0d0 + fy3(:) = 0.0d0 + fz3(:) = 0.0d0 +!----------------------------------------------------------------------------------------- +! triple loop for the 2 and 3-body forces +! do 20 i +! do 30 j +! do 40 k +! the pairlists lsta and lstb are used for the perodic boundary conditions +!----------------------------------------------------------------------------------------- + + DO i = 1, nat + CALL sw_subfeniat_l(i, nat, nnbrx, rel, p, p3, fx, fy, fz, fx3, fy3, fz3, lstb, lsta, isigma, sigma) + END DO +!----------------------------------------------------------------------------------------- +!*****if necessary,******** + DO i = 1, nat + fx(i) = fx(i) + fx3(i) + fy(i) = fy(i) + fy3(i) + fz(i) = fz(i) + fz3(i) + END DO +!----------------------------------------------------------------------------------------- + + DO i = 1, nat + fxyz(1, i) = fx(i)*esigma + fxyz(2, i) = fy(i)*esigma + fxyz(3, i) = fz(i)*esigma + END DO + etot = (p + p3)*eps + DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel) + END SUBROUTINE eip_stillinger_weber_silicon +!End of the force calculation------------------------------------------------------------- + +! ************************************************************************************************** +!> \brief ... +!> \param c ... +!> \return ... +! ************************************************************************************************** + REAL(8) FUNCTION f(c) + REAL(8) :: c + + REAL(8) :: aa, bb, c4, crainv, ra + + PARAMETER(aa=7.049556277d0) + PARAMETER(ra=1.8d0, bb=0.6022245584d0) + + IF ((c - ra) < 0.d0) THEN + crainv = 1.0d0/(c - ra) + c4 = c*c*c*c + f = aa*bb*4.0d0/(c4*c)*dexp(crainv) + aa*(bb/(c4) - 1.0d0)*dexp(crainv)*crainv*crainv + ELSE + f = 0.d0 + END IF + + RETURN + END FUNCTION f + +! ************************************************************************************************** +!> \brief ... +!> \param d ... +!> \return ... +! ************************************************************************************************** + REAL(8) FUNCTION pe(d) + REAL(8) :: d + + REAL(8), PARAMETER :: aa = 7.049556277d0, bb = 0.6022245584d0, & + ra = 1.8d0 + + IF ((d - ra) < 0.d0) THEN + pe = aa*(bb/(d*d*d*d) - 1.0d0)*dexp(1.0d0/(d - ra)) + ELSE + pe = 0.d0 + END IF + RETURN + END FUNCTION pe + +!------------------------------------------------------------------------------------------ +! ************************************************************************************************** +!> \brief ... +!> \param iat ... +!> \param nn ... +!> \param ncx ... +!> \param ll1 ... +!> \param ll2 ... +!> \param ll3 ... +!> \param l1 ... +!> \param l2 ... +!> \param l3 ... +!> \param myspace ... +!> \param rxyz ... +!> \param icell ... +!> \param lstb ... +!> \param lay ... +!> \param rel ... +!> \param cut2 ... +!> \param indlst ... +! ************************************************************************************************** + SUBROUTINE sw_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, & + rxyz, icell, lstb, lay, rel, cut2, indlst) +! finds the neighbours of atom iat (specified by lsta and lstb) and and +! the relative position rel of iat with respect to these neighbours + INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, & + myspace + REAL(8) :: rxyz(3, nn) + INTEGER :: icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), lstb(0:myspace - 1), lay(nn) + REAL(8) :: rel(5, 0:myspace - 1), cut2 + INTEGER :: indlst + + INTEGER :: jat, jj, k1, k2, k3 + REAL(8) :: rr2, tt, tti, xrel, yrel, zrel + + DO k3 = l3 - 1, l3 + 1 + DO k2 = l2 - 1, l2 + 1 + DO k1 = l1 - 1, l1 + 1 + DO jj = 1, icell(0, k1, k2, k3) + jat = icell(jj, k1, k2, k3) + IF (jat == iat) CYCLE + xrel = rxyz(1, iat) - rxyz(1, jat) + yrel = rxyz(2, iat) - rxyz(2, jat) + zrel = rxyz(3, iat) - rxyz(3, jat) + rr2 = xrel**2 + yrel**2 + zrel**2 + IF (rr2 <= cut2) THEN + indlst = MIN(indlst, myspace - 1) + lstb(indlst) = lay(jat) +! write(6,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat) + tt = SQRT(rr2) + tti = 1.d0/tt + rel(1, indlst) = xrel*tti + rel(2, indlst) = yrel*tti + rel(3, indlst) = zrel*tti + rel(4, indlst) = tt + rel(5, indlst) = tti + indlst = indlst + 1 + END IF + END DO + END DO + END DO + END DO + + RETURN + END SUBROUTINE sw_sublstiat_l + +! ************************************************************************************************** +!> \brief ... +!> \param i ... +!> \param nat ... +!> \param nnbrx ... +!> \param rel ... +!> \param p ... +!> \param p3 ... +!> \param fx ... +!> \param fy ... +!> \param fz ... +!> \param fx3 ... +!> \param fy3 ... +!> \param fz3 ... +!> \param lstb ... +!> \param lsta ... +!> \param isigma ... +!> \param sigma ... +! ************************************************************************************************** + SUBROUTINE sw_subfeniat_l(i, nat, nnbrx, rel, p, p3, fx, fy, fz, fx3, fy3, fz3, lstb, lsta, isigma, sigma) + INTEGER, INTENT(IN) :: i, nat, nnbrx + REAL(8), INTENT(IN) :: rel(5, nnbrx*nat) + REAL(8), INTENT(INOUT) :: p, p3, fx(nat), fy(nat), fz(nat), & + fx3(nat), fy3(nat), fz3(nat) + INTEGER, INTENT(IN) :: lstb(nnbrx*nat), lsta(2, nat) + REAL(8), INTENT(IN) :: isigma, sigma + + INTEGER :: Ipb, Ipe, j, k, l, m, nij + REAL(8) :: aa, bb, c4, cosijk, cosijk3, cosikj, cosikj3, cosjik, cosjik3, crainv, force, & + gam, hi, hixij, hixij0, hixij1, hixik, hixik0, hixik1, hiyij, hiyij0, hiyij1, hiyik, & + hiyik0, hiyik1, HIZIJ, HIZIJ0, HIZIJ1, HIZIK, HIZIK0, HIZIK1, hj, hjxij, hjxij0, hjxij1, & + hjxjk, hjxjk0, hjxjk1, hjyij, hjyij0, hjyij1, hjyjk, hjyjk0, hjyjk1, HJZIJ, HJZIJ0, & + HJZIJ1, HJZJK, HJZJK0, HJZJK1, hk, hkxik, hkxik0, HKXIK1, hkxkj, hkxkj0, hkxkj1, hkyik, & + hkyik0, hkyik1, hkykj, hkykj0, hkykj1, HKZIK, HKZIK0, HKZIK1, HKZKJ, HKZKJ0, HKZKJ1, & + invrij, invrija, invrik, invrika, invrjk, invrjka, ra, ramda, refi, refj + REAL(8) :: refk, rij, rija, rik, rika, rjk, rjka, xij, xik, xjk, yij, yik, yjk, zij, zik, zjk + + PARAMETER(gam=1.2d0, ramda=21.0d0, aa=7.049556277d0) + PARAMETER(ra=1.8d0, bb=0.6022245584d0) + + Ipb = lsta(1, i) + Ipe = lsta(2, i) + + DO l = Ipb, Ipe + j = lstb(l) + IF (j <= i) CYCLE + nij = 0 + rij = rel(4, l)*isigma + invrij = rel(5, l)*sigma + xij = rel(1, l)*rij + yij = rel(2, l)*rij + zij = rel(3, l)*rij + + IF (rij >= 2.d0*ra) CYCLE + IF (rij < ra) THEN + crainv = 1.0d0/(rij - ra) + c4 = rij*rij*rij*rij + force = aa*bb*4.0d0/(c4*rij)*dexp(crainv) + aa*(bb/(c4) - 1.0d0)*dexp(crainv)*crainv*crainv + + fx(i) = force*xij*invrij + fx(i) + fy(i) = force*yij*invrij + fy(i) + fz(i) = force*zij*invrij + fz(i) + + fx(j) = -force*xij*invrij + fx(j) + fy(j) = -force*yij*invrij + fy(j) + fz(j) = -force*zij*invrij + fz(j) + + p = p + aa*(bb/(rij*rij*rij*rij) - 1.0d0)*dexp(1.0d0/(rij - ra)) + + nij = 1 + END IF + + DO m = Ipb, Ipe + k = lstb(m) + IF (k <= j) CYCLE + invrik = rel(5, m)*sigma + rik = rel(4, m)*isigma + xik = rel(1, m)*rik + yik = rel(2, m)*rik + zik = rel(3, m)*rik + + IF ((rik >= ra) .AND. (nij == 0)) CYCLE + + xjk = xik - xij + yjk = yik - yij + zjk = zik - zij + + rjk = SQRT(xjk*xjk + yjk*yjk + zjk*zjk) + invrjk = 1.d0/rjk + + IF ((rjk >= ra) .AND. (nij == 0)) CYCLE + cosjik = (xij*xik + yij*yik + zij*zik)*(invrij*invrik) + cosijk = (-xij*xjk - yij*yjk - zij*zjk)*(invrij*invrjk) + cosikj = (xik*xjk + yik*yjk + zik*zjk)*(invrik*invrjk) + cosjik3 = cosjik + 1.0d0/3.0d0 + cosijk3 = cosijk + 1.0d0/3.0d0 + cosikj3 = cosikj + 1.0d0/3.0d0 + + rija = rij - ra + rika = rik - ra + rjka = rjk - ra + + invrija = 1.d0/rija + invrika = 1.d0/rika + invrjka = 1.d0/rjka + + IF (rija >= 0.0d0) THEN + refi = 0.0d0 + refj = 0.0d0 + refk = ramda*EXP(gam*invrika + gam*invrjka) + ELSE IF ((rija < 0.0d0) .AND. (rika < 0.0d0)) THEN + IF (rjka < 0.0d0) THEN + refi = ramda*EXP(gam*invrija + gam*invrika) + refj = ramda*EXP(gam*invrija + gam*invrjka) + refk = ramda*EXP(gam*invrika + gam*invrjka) + ELSE + refi = ramda*EXP(gam*invrija + gam*invrika) + refj = 0.0d0 + refk = 0.0d0 + END IF + ELSE IF ((rija < 0.0d0) .AND. (rjka < 0.0d0)) THEN + refi = 0.0d0 + refj = ramda*EXP(gam*invrija + gam*invrjka) + refk = 0.0d0 + ELSE + CYCLE + END IF + + hi = refi*cosjik3*cosjik3 + hj = refj*cosijk3*cosijk3 + hk = refk*cosikj3*cosikj3 + p3 = p3 + hi + hj + hk + + hixij0 = 2.0d0*(xik*invrik - xij*cosjik*invrij) + hixij1 = gam*xij*cosjik3*(invrija*invrija) + hixij = refi*cosjik3*(hixij0 - hixij1)*invrij + hixik0 = 2.0d0*(xij*invrij - xik*cosjik*invrik) + hixik1 = gam*xik*cosjik3*(invrika*invrika) + hixik = refi*cosjik3*(hixik0 - hixik1)*invrik + hjxij0 = 2.0d0*(-xjk*invrjk - xij*cosijk*invrij) + hjxij1 = gam*xij*cosijk3*(invrija*invrija) + hjxij = refj*cosijk3*(hjxij0 - hjxij1)*invrij + hkxik0 = 2.0d0*(xjk*invrjk - xik*cosikj*invrik) + hkxik1 = gam*xik*cosikj3*(invrika*invrika) + hkxik = refk*cosikj3*(hkxik0 - hkxik1)*invrik + hjxjk0 = 2.0d0*(-xij*invrij - xjk*cosijk*invrjk) + hjxjk1 = gam*xjk*cosijk3*(invrjka*invrjka) + hjxjk = refj*cosijk3*(hjxjk0 - hjxjk1)*invrjk + hkxkj0 = 2.0d0*(-xik*invrik + xjk*cosikj*invrjk) + hkxkj1 = gam*xjk*cosikj3*(invrjka*invrjka) + hkxkj = refk*cosikj3*(hkxkj0 + hkxkj1)*invrjk + + hiyij0 = 2.0d0*(yik*invrik - yij*cosjik*invrij) + hiyij1 = gam*yij*cosjik3*(invrija*invrija) + hiyij = refi*cosjik3*(hiyij0 - hiyij1)*invrij + hiyik0 = 2.0d0*(yij*invrij - yik*cosjik*invrik) + hiyik1 = gam*yik*cosjik3*(invrika*invrika) + hiyik = refi*cosjik3*(hiyik0 - hiyik1)*invrik + hjyij0 = 2.0d0*(-yjk*invrjk - yij*cosijk*invrij) + hjyij1 = gam*yij*cosijk3*(invrija*invrija) + hjyij = refj*cosijk3*(hjyij0 - hjyij1)*invrij + hkyik0 = 2.0d0*(yjk*invrjk - yik*cosikj*invrik) + hkyik1 = gam*yik*cosikj3*(invrika*invrika) + hkyik = refk*cosikj3*(hkyik0 - hkyik1)*invrik + hjyjk0 = 2.0d0*(-yij*invrij - yjk*cosijk*invrjk) + hjyjk1 = gam*yjk*cosijk3*(invrjka*invrjka) + hjyjk = refj*cosijk3*(hjyjk0 - hjyjk1)*invrjk + hkykj0 = 2.0d0*(-yik*invrik + yjk*cosikj*invrjk) + hkykj1 = gam*yjk*cosikj3*(invrjka*invrjka) + hkykj = refk*cosikj3*(hkykj0 + hkykj1)*invrjk + + hizij0 = 2.0d0*(zik*invrik - zij*cosjik*invrij) + hizij1 = gam*zij*cosjik3*(invrija*invrija) + hizij = refi*cosjik3*(hizij0 - hizij1)*invrij + hizik0 = 2.0d0*(zij*invrij - zik*cosjik*invrik) + hizik1 = gam*zik*cosjik3*(invrika*invrika) + hizik = refi*cosjik3*(hizik0 - hizik1)*invrik + hjzij0 = 2.0d0*(-zjk*invrjk - zij*cosijk*invrij) + hjzij1 = gam*zij*cosijk3*(invrija*invrija) + hjzij = refj*cosijk3*(hjzij0 - hjzij1)*invrij + hkzik0 = 2.0d0*(zjk*invrjk - zik*cosikj*invrik) + hkzik1 = gam*zik*cosikj3*(invrika*invrika) + hkzik = refk*cosikj3*(hkzik0 - hkzik1)*invrik + hjzjk0 = 2.0d0*(-zij*invrij - zjk*cosijk*invrjk) + hjzjk1 = gam*zjk*cosijk3*(invrjka*invrjka) + hjzjk = refj*cosijk3*(hjzjk0 - hjzjk1)*invrjk + hkzkj0 = 2.0d0*(-zik*invrik + zjk*cosikj*invrjk) + hkzkj1 = gam*zjk*cosikj3*(invrjka*invrjka) + hkzkj = refk*cosikj3*(hkzkj0 + hkzkj1)*invrjk + + fx3(i) = fx3(i) - hixij - hixik - hjxij - hkxik + fy3(i) = fy3(i) - hiyij - hiyik - hjyij - hkyik + fz3(i) = fz3(i) - hizij - hizik - hjzij - hkzik + + fx3(j) = fx3(j) + hixij + hjxij - hjxjk + hkxkj + fy3(j) = fy3(j) + hiyij + hjyij - hjyjk + hkykj + fz3(j) = fz3(j) + hizij + hjzij - hjzjk + hkzkj + + fx3(k) = fx3(k) + hixik + hkxik - hkxkj + hjxjk + fy3(k) = fy3(k) + hiyik + hkyik - hkykj + hjyjk + fz3(k) = fz3(k) + hizik + hkzik - hkzkj + hjzjk + END DO + END DO + END SUBROUTINE sw_subfeniat_l + +! ************************************************************************************************** +!> \brief ... +!> \param nat ... +!> \param alat ... +!> \param rxyz ... +!> \param fxyz ... +!> \param etot ... +!> \param count ... +! ************************************************************************************************** + SUBROUTINE eip_tersoff_silicon(nat, alat, rxyz, fxyz, etot, count) +!***************************************************************************************** +! This subroutine evaluates the Tersoff Silicon potential with linear scaling +! COPYRIGHT +! Copyright (C) 2009 AIST, UNIBAS +! This file is distributed under the terms of the +! GNU General Public License, see +! http://www.gnu.org/copyleft/gpl.txt . +! +! Implementation: Original version was written by Kengo Nishio, AIST Tsukuba (JP) +! Improved by M. Amsler, S. Goedecker, Basel University (CH), 2009 +! +! Note: +! Parameters and functional form from PRL 61, 2879 (1988) and PRB 39, 5566 (1989) +! +! Input: +! nat, integer: the number of atoms +! alat, real(8), dim(3) : the three edges of the orthoromic simulation cell, periodic boundaries are applied +! and atoms outside the cell will be brought back into the box +! rxyz, real(8), dim(3,nat) : the xyz cartesian components of the atomic positions in Angstroem +! +! Output: +! fxyz, real(8), dim(3,nat): the xyz cartesian forces in on the corresponding atomic components n eV/A +! etot, real(8) : total potential energy, 2-body and 3-body, in eV +! count, real(8) : increased by 1.d0 at each call of this subroutine, +! needs to be initialized to 0.d0 before calling this routine for the first time +!***************************************************************************************** + INTEGER :: nat + REAL(8) :: alat(3), rxyz(3, nat), fxyz(3, nat), & + etot, count + + INTEGER :: iat, NNmax, Npmax + INTEGER, ALLOCATABLE, DIMENSION(:) :: Kinds, lstb + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta + REAL(8) :: Uatot, Urtot, xbox, ybox, zbox + REAL(8), ALLOCATABLE, DIMENSION(:) :: dkEij, UadUrdf, XYZRrefdf + REAL(8), DIMENSION(1:2) :: bcsq, Co_bcd, dsq, h, Pmass, Pn + REAL(8), DIMENSION(1:2, 1:2) :: ala, alr, Ca, Cr, R1, R2, X + + INTEGER:: nnbrx, nnbrxt + INTEGER :: i + + count = count + 1.d0 + + DO iat = 1, nat + rxyz(1, iat) = MODULO(MODULO(rxyz(1, iat), alat(1)), alat(1)) + rxyz(2, iat) = MODULO(MODULO(rxyz(2, iat), alat(2)), alat(2)) + rxyz(3, iat) = MODULO(MODULO(rxyz(3, iat), alat(3)), alat(3)) + END DO + + ALLOCATE (Kinds(1:nat)) + + nnbrx = 24 + nnbrxt = 3*nnbrx/2 + nnmax = nnbrxt*nat + npmax = nnbrxt*nat + ALLOCATE (lsta(2, nat), lstb(nnbrxt*nat)) + ALLOCATE (XYZRrefdf(1:6*Npmax), UadUrdf(1:3*Npmax), dkEij(1:3*NNmax)) + + DO i = 1, nat + kinds(i) = 2 !Since all atoms are Si, all of kind 2 + END DO + fxyz = 0.0d0 + xbox = alat(1); ybox = alat(2); zbox = alat(3) + CALL tersoff_parameters(R1, R2, Cr, Ca, alr, ala, X, Pn, Co_bcd, bcsq, dsq, h, Pmass) + CALL tersoff_pairlist_energy_forces(nat, Npmax, NNmax, xbox, ybox, zbox, Kinds, rxyz, R1, R2, Cr, & + Ca, alr, ala, X, XYZRrefdf, UadUrdf, Urtot, lsta, lstb, nnbrx, & + Pn, Co_bcd, bcsq, dsq, h, fxyz, Uatot, dkEij) + etot = Urtot + Uatot + DEALLOCATE (Kinds, XYZRrefdf, UadUrdf, dkEij, lsta, lstb) + END SUBROUTINE eip_tersoff_silicon +!----------------------------------------------------------------------------------------- +! ************************************************************************************************** +!> \brief ... +!> \param R1 ... +!> \param R2 ... +!> \param Cr ... +!> \param Ca ... +!> \param alr ... +!> \param ala ... +!> \param X ... +!> \param Pn ... +!> \param Co_bcd ... +!> \param bcsq ... +!> \param dsq ... +!> \param h ... +!> \param Pmass ... +! ************************************************************************************************** + SUBROUTINE tersoff_parameters(R1, R2, Cr, Ca, alr, ala, X, Pn, Co_bcd, bcsq, dsq, h, Pmass) + + REAL(8), DIMENSION(1:2, 1:2), INTENT(out) :: R1, R2, Cr, Ca, alr, ala, X + REAL(8), DIMENSION(1:2), INTENT(out) :: Pn, Co_bcd, bcsq, dsq, h, Pmass + + REAL(8), PARAMETER :: C_ala = 2.2119d0, C_alr = 3.4879d0, C_b = 1.5724d-7, C_c = 3.8049d4, & + C_Ca = 3.4674d2, C_Cr = 1.3936d3, C_d = 4.3484d0, C_h = -5.7058d-1, C_mass = 12.0d0, & + C_n = 7.2751d-1, C_R1 = 1.8d0, C_R2 = 2.1d0, Si_ala = 1.7322d0, Si_alr = 2.4799d0, & + Si_b = 1.1000d-6, Si_c = 1.0039d5, Si_Ca = 4.7118d2, Si_Cr = 1.8308d3, Si_d = 1.6217d1, & + Si_h = -5.9825d-1, Si_mass = 28.0855d0, Si_n = 7.8734d-1, Si_R1 = 2.7d0, Si_R2 = 3.3d0 + +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Parameter for carbon, not used in this version +!Increased Cutoff, originally 3.0d0 + + Cr(1, 1) = C_Cr + Cr(2, 2) = Si_Cr + Cr(1, 2) = dsqrt(Cr(1, 1)*Cr(2, 2)) + Cr(2, 1) = Cr(1, 2) + + Ca(1, 1) = C_Ca + Ca(2, 2) = Si_Ca + Ca(1, 2) = dsqrt(Ca(1, 1)*Ca(2, 2)) + Ca(2, 1) = Ca(1, 2) + + R1(1, 1) = C_R1 + R1(2, 2) = Si_R1 + R1(1, 2) = dsqrt(R1(1, 1)*R1(2, 2)) + R1(2, 1) = R1(1, 2) + + R2(1, 1) = C_R2 + R2(2, 2) = Si_R2 + R2(1, 2) = dsqrt(R2(1, 1)*R2(2, 2)) + R2(2, 1) = R2(1, 2) + + X(1, 1) = 1.0d0 + X(2, 2) = 1.0d0 + X(1, 2) = 0.9776d0 + X(2, 1) = 0.9776d0 + + alr(1, 1) = C_alr + alr(2, 2) = Si_alr + alr(1, 2) = 0.5d0*(alr(1, 1) + alr(2, 2)) + alr(2, 1) = alr(1, 2) + + ala(1, 1) = C_ala + ala(2, 2) = Si_ala + ala(1, 2) = 0.5d0*(ala(1, 1) + ala(2, 2)) + ala(2, 1) = ala(1, 2) + + Pn(1) = C_n + Pn(2) = Si_n + + Co_bcd(1) = C_b*(1.0d0 + C_c*C_c/(C_d*C_d)) + Co_bcd(2) = Si_b*(1.0d0 + Si_c*Si_c/(Si_d*Si_d)) + + bcsq(1) = C_b*C_c*C_c + bcsq(2) = Si_b*Si_c*Si_c + + dsq(1) = C_d*C_d + dsq(2) = Si_d*Si_d + + h(1) = C_h + h(2) = Si_h + + Pmass(1) = C_mass + Pmass(2) = Si_mass + + RETURN + END SUBROUTINE tersoff_parameters +!----------------------------------------------------------------------------------------- +! ************************************************************************************************** +!> \brief ... +!> \param Nmol ... +!> \param Npmax ... +!> \param NNmax ... +!> \param xbox ... +!> \param ybox ... +!> \param zbox ... +!> \param Kinds ... +!> \param R ... +!> \param R1 ... +!> \param R2 ... +!> \param Cr ... +!> \param Ca ... +!> \param alr ... +!> \param ala ... +!> \param X ... +!> \param XYZRrefdf ... +!> \param UadUrdf ... +!> \param Urtot ... +!> \param lsta ... +!> \param lstb ... +!> \param nnbrx ... +!> \param Pn ... +!> \param Co_bcd ... +!> \param bcsq ... +!> \param dsq ... +!> \param h ... +!> \param F ... +!> \param Uatot ... +!> \param dkEij ... +! ************************************************************************************************** + SUBROUTINE tersoff_pairlist_energy_forces(Nmol, Npmax, NNmax, xbox, ybox, zbox, Kinds, R, R1, R2, Cr, Ca, alr, ala, X, & + XYZRrefdf, UadUrdf, Urtot, lsta, lstb, nnbrx, & + Pn, Co_bcd, bcsq, dsq, h, F, Uatot, dkEij) + INTEGER, INTENT(in) :: Nmol, Npmax, NNmax + REAL(8), INTENT(in) :: xbox, ybox, zbox + INTEGER, DIMENSION(1:Nmol), INTENT(in) :: Kinds + REAL(8), DIMENSION(1:3*Nmol), INTENT(in) :: R + REAL(8), DIMENSION(1:2, 1:2), INTENT(in) :: R1, R2, Cr, Ca, alr, ala, X + REAL(8), DIMENSION(1:6*Npmax), INTENT(out) :: XYZRrefdf + REAL(8), DIMENSION(1:3*Npmax), INTENT(out) :: UadUrdf + REAL(8), INTENT(out) :: Urtot + INTEGER :: lsta(2, Nmol) + INTEGER, INTENT(inout) :: nnbrx + INTEGER :: lstb(nnbrx*Nmol) + REAL(8), DIMENSION(1:2), INTENT(in) :: Pn, Co_bcd, bcsq, dsq, h + REAL(8), DIMENSION(1:3*Nmol), INTENT(out) :: F + REAL(8), INTENT(out) :: Uatot + REAL(8), DIMENSION(1:3*NNmax) :: dkEij + + INTEGER :: i, iam, iat, ii, il, in, indlst, indlstx, Ipb, istopg, jat, l1, l2, l3, laymx, & + ll1, ll2, ll3, myspace, myspaceout, nat, ncx, ndat, nn, npjkx, npjx, npr, Nptot + INTEGER, ALLOCATABLE, DIMENSION(:) :: lay + INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell + REAL(8) :: alat(3), cut, cut2, pi, rlc1i, rlc2i, & + rlc3i, rxyz0(3, Nmol), xhalf, yhalf, & + zhalf + REAL(8), ALLOCATABLE, DIMENSION(:, :) :: rel, rxyz + + pi = dacos(-1.0d0) + + xhalf = 0.5d0*xbox + yhalf = 0.5d0*ybox + zhalf = 0.5d0*zbox + + nat = Nmol + + alat(1) = xbox + alat(2) = ybox + alat(3) = zbox + + DO iat = 1, nat + jat = 3*(iat - 1) + rxyz0(1, iat) = R(jat + 1) + rxyz0(2, iat) = R(jat + 2) + rxyz0(3, iat) = R(jat + 3) + END DO + + cut = R2(2, 2) - 1.d-9 + +! linear scaling calculation of verlet list + ll1 = INT(alat(1)/cut) + IF (ll1 < 1) CPABORT("alat(1) too small") + ll2 = INT(alat(2)/cut) + IF (ll2 < 1) CPABORT("alat(2) too small") + ll3 = INT(alat(3)/cut) + IF (ll3 < 1) CPABORT("alat(3) too small") + +! determine number of threadsi (this version is only singlethreaded) + npr = 1 +! linear scaling calculation of verlet list + + ncx = 8 + DO + ncx = ncx*2 + ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3)) + icell(0, :, :, :) = 0 + rlc1i = ll1/alat(1) + rlc2i = ll2/alat(2) + rlc3i = ll3/alat(3) + + DO iat = 1, nat + l1 = INT(rxyz0(1, iat)*rlc1i) + l2 = INT(rxyz0(2, iat)*rlc2i) + l3 = INT(rxyz0(3, iat)*rlc3i) + + ii = icell(0, l1, l2, l3) + ii = ii + 1 + icell(0, l1, l2, l3) = ii + IF (ii > ncx) THEN + DEALLOCATE (icell) + EXIT + END IF + icell(ii, l1, l2, l3) = iat + END DO + IF (ALLOCATED(icell)) EXIT + END DO + +! duplicate all atoms within boundary layer + laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8) + nn = nat + laymx + ALLOCATE (rxyz(3, nn), lay(nn)) + DO iat = 1, nat + lay(iat) = iat + rxyz(1, iat) = rxyz0(1, iat) + rxyz(2, iat) = rxyz0(2, iat) + rxyz(3, iat) = rxyz0(3, iat) + END DO + il = nat +! xy plane + DO l2 = 0, ll2 - 1 + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, l2, 0) + icell(0, l1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, l1, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, l2, ll3 - 1) + icell(0, l1, l2, -1) = in + DO ii = 1, in + i = icell(ii, l1, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + END DO + +! yz plane + DO l3 = 0, ll3 - 1 + DO l2 = 0, ll2 - 1 + + in = icell(0, 0, l2, l3) + icell(0, ll1, l2, l3) = in + DO ii = 1, in + i = icell(ii, 0, l2, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, l2, l3) + icell(0, -1, l2, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + END DO + +! xz plane + DO l3 = 0, ll3 - 1 + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, 0, l3) + icell(0, l1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, l1, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, l1, ll2 - 1, l3) + icell(0, l1, -1, l3) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + END DO + +! x axis + DO l1 = 0, ll1 - 1 + + in = icell(0, l1, 0, 0) + icell(0, l1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, l1, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, 0, ll3 - 1) + icell(0, l1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, l1, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, l1, ll2 - 1, 0) + icell(0, l1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, l1, ll2 - 1, ll3 - 1) + icell(0, l1, -1, -1) = in + DO ii = 1, in + i = icell(ii, l1, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, l1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + +! y axis + DO l2 = 0, ll2 - 1 + + in = icell(0, 0, l2, 0) + icell(0, ll1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, 0, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, l2, ll3 - 1) + icell(0, ll1, l2, -1) = in + DO ii = 1, in + i = icell(ii, 0, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, l2, 0) + icell(0, -1, l2, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, l2, ll3 - 1) + icell(0, -1, l2, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, l2, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, l2, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + END DO + +! z axis + DO l3 = 0, ll3 - 1 + + in = icell(0, 0, 0, l3) + icell(0, ll1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, 0, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, 0, l3) + icell(0, -1, ll2, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, 0, ll2 - 1, l3) + icell(0, ll1, -1, l3) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, l3) + icell(0, -1, -1, l3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, l3) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, l3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + END DO + + END DO + +! corners + in = icell(0, 0, 0, 0) + icell(0, ll1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, 0, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, 0, 0) + icell(0, -1, ll2, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, ll2 - 1, 0) + icell(0, ll1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, 0) + icell(0, -1, -1, ll3) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, 0) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, ll3) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) + alat(3) + END DO + + in = icell(0, 0, 0, ll3 - 1) + icell(0, ll1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, 0, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, 0, ll3 - 1) + icell(0, -1, ll2, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, 0, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, ll2, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) + alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, 0, ll2 - 1, ll3 - 1) + icell(0, ll1, -1, -1) = in + DO ii = 1, in + i = icell(ii, 0, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, ll1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) + alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1) + icell(0, -1, -1, -1) = in + DO ii = 1, in + i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1) + il = il + 1 + IF (il > nn) CPABORT("enlarge laymx") + lay(il) = i + icell(ii, -1, -1, -1) = il + rxyz(1, il) = rxyz(1, i) - alat(1) + rxyz(2, il) = rxyz(2, i) - alat(2) + rxyz(3, il) = rxyz(3, i) - alat(3) + END DO + + nnbrx = 3*nnbrx/2 + ALLOCATE (rel(5, nnbrx*nat)) + indlstx = 0 + + npr = 1 + iam = 0 + + cut2 = cut**2 +! assign contiguous portions of the arrays lstb and rel to the threads + myspace = (nat*nnbrx)/npr + IF (iam == 0) myspaceout = myspace +! Verlet list, relative positions + indlst = 0 + DO l3 = 0, ll3 - 1 + DO l2 = 0, ll2 - 1 + DO l1 = 0, ll1 - 1 + DO ii = 1, icell(0, l1, l2, l3) + iat = icell(ii, l1, l2, l3) + IF (((iat - 1)*npr)/nat == iam) THEN +! write(6,*) 'sublstiat:iam,iat',iam,iat + lsta(1, iat) = iam*myspace + indlst + 1 + CALL tersoff_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, & + rxyz, icell, lstb(iam*myspace + 1), lay, rel(1, iam*myspace + 1), cut2, indlst) + lsta(2, iat) = iam*myspace + indlst + ipb = lsta(1, iat) + ndat = lsta(2, iat) - lsta(1, iat) + 1 + END IF + + END DO + END DO + END DO + END DO + indlstx = MAX(indlstx, indlst) + + IF (indlstx >= myspaceout) CPABORT("NNBRX too small") + npr = 1 + iam = 0 + + npjx = 300; npjkx = 6000 + istopg = 0 +!end of creating pairlist part------------------------------------------------------------ +!Energy----------------------------------------------------------------------------------- + Urtot = 0.0d0 + Nptot = 0 + + F = 0.0d0 + Uatot = 0.0d0 + DO_I: DO i = 1, Nmol + CALL tersoff_subeniat_l(i, Nmol, Npmax, Kinds, X, R1, R2, Cr, Ca, alr, ala, XYZRrefdf, UadUrdf, Urtot, lsta, lstb, nnbrx, rel, pi) + END DO DO_I + + Urtot = 0.5d0*Urtot +!Force------------------------------------------------------------------------------------ + F = 0.0d0 + Uatot = 0.0d0 + + DO_If: DO i = 1, Nmol + CALL tersoff_subfiat_l(i,Nmol,Npmax,NNmax,Kinds,Pn,Co_bcd,bcsq,dsq,h,XYZRrefdf,UadUrdf,F,Uatot,dkEij,lsta,lstb,nnbrx) + END DO DO_If + + F = 0.5d0*F + Uatot = 0.5d0*Uatot +!----------------------------------------------------------------------------------------- + DEALLOCATE (rxyz, icell, lay, rel) + RETURN + END SUBROUTINE tersoff_pairlist_energy_forces + +!----------------------------------------------------------------------------------------- +! ************************************************************************************************** +!> \brief ... +!> \param iat ... +!> \param nn ... +!> \param ncx ... +!> \param ll1 ... +!> \param ll2 ... +!> \param ll3 ... +!> \param l1 ... +!> \param l2 ... +!> \param l3 ... +!> \param myspace ... +!> \param rxyz ... +!> \param icell ... +!> \param lstb ... +!> \param lay ... +!> \param rel ... +!> \param cut2 ... +!> \param indlst ... +! ************************************************************************************************** + SUBROUTINE tersoff_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, & + rxyz, icell, lstb, lay, rel, cut2, indlst) +! finds the neighbours of atom iat (specified by lsta and lstb) and and +! the relative position rel of iat with respect to these neighbours + INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, & + myspace + REAL(8) :: rxyz(3, nn) + INTEGER :: icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), lstb(0:myspace - 1), lay(nn) + REAL(8) :: rel(5, 0:myspace - 1), cut2 + INTEGER :: indlst + + INTEGER :: jat, jj, k1, k2, k3 + REAL(8) :: rr2, tt, tti, xrel, yrel, zrel + + DO k3 = l3 - 1, l3 + 1 + DO k2 = l2 - 1, l2 + 1 + DO k1 = l1 - 1, l1 + 1 + DO jj = 1, icell(0, k1, k2, k3) + jat = icell(jj, k1, k2, k3) + IF (jat == iat) CYCLE + xrel = rxyz(1, iat) - rxyz(1, jat) + yrel = rxyz(2, iat) - rxyz(2, jat) + zrel = rxyz(3, iat) - rxyz(3, jat) + rr2 = xrel**2 + yrel**2 + zrel**2 + IF (rr2 <= cut2) THEN + indlst = MIN(indlst, myspace - 1) + lstb(indlst) = lay(jat) +! write(6,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat) + tt = SQRT(rr2) + tti = 1.d0/tt + rel(1, indlst) = xrel*tti + rel(2, indlst) = yrel*tti + rel(3, indlst) = zrel*tti + rel(4, indlst) = tt + rel(5, indlst) = tti + indlst = indlst + 1 + END IF + END DO + END DO + END DO + END DO + + RETURN + END SUBROUTINE tersoff_sublstiat_l + +! ************************************************************************************************** +!> \brief ... +!> \param i ... +!> \param Nmol ... +!> \param Npmax ... +!> \param Kinds ... +!> \param X ... +!> \param R1 ... +!> \param R2 ... +!> \param Cr ... +!> \param Ca ... +!> \param alr ... +!> \param ala ... +!> \param XYZRrefdf ... +!> \param UadUrdf ... +!> \param Urtot ... +!> \param lsta ... +!> \param lstb ... +!> \param nnbrx ... +!> \param rel ... +!> \param pi ... +! ************************************************************************************************** + SUBROUTINE tersoff_subeniat_l(i,Nmol,Npmax,Kinds,X,R1,R2,Cr,Ca,alr,ala,XYZRrefdf,UadUrdf,Urtot,lsta,lstb,nnbrx,rel,pi) + INTEGER :: i + INTEGER, INTENT(in) :: Nmol, Npmax + INTEGER, DIMENSION(1:Nmol), INTENT(in) :: Kinds + REAL(8), DIMENSION(1:2, 1:2), INTENT(in) :: X, R1, R2, Cr, Ca, alr, ala + REAL(8), DIMENSION(1:6*Npmax), INTENT(inout) :: XYZRrefdf + REAL(8), DIMENSION(1:3*Npmax), INTENT(inout) :: UadUrdf + REAL(8), INTENT(inout) :: Urtot + INTEGER, INTENT(in) :: lsta(2, Nmol), nnbrx, lstb(nnbrx*Nmol) + REAL(8), INTENT(in) :: rel(5, nnbrx*Nmol) + REAL(8) :: pi + + INTEGER :: j, Ki, Kj, l, Nppt3, Nppt6, Nptot + REAL(8) :: alaij, alrij, dfij, fij, PL1, PL2, R1ij, & + R2ij, Rij, Rreij, Ua, Ur, Xij, Yij, Zij + +! ####################################### +! # Calculate XYZRrefdf, UadUrdf, Urtot # +! ####################################### + Ki = Kinds(i) + + DO_J: DO l = lsta(1, i), lsta(2, i) + j = lstb(l) + + Kj = Kinds(j) + R2ij = R2(Ki, Kj) + Rij = rel(4, l) + Xij = rel(1, l) + Yij = rel(2, l) + Zij = rel(3, l) + Nptot = l + + Nppt3 = 3*(Nptot - 1) + Nppt6 = 6*(Nptot - 1) + Rreij = rel(5, l) + + XYZRrefdf(Nppt6 + 1) = Xij + XYZRrefdf(Nppt6 + 2) = Yij + XYZRrefdf(Nppt6 + 3) = Zij + XYZRrefdf(Nppt6 + 4) = Rreij + + alrij = alr(Ki, Kj) + alaij = ala(Ki, Kj) + + Ur = Cr(Ki, Kj)*dexp(-alrij*Rij) + Ua = -Ca(Ki, Kj)*dexp(-alaij*Rij)*X(Ki, Kj) + R1ij = R1(Ki, Kj) + + IF (Rij <= R1ij) THEN + XYZRrefdf(Nppt6 + 5) = 1.0d0 + XYZRrefdf(Nppt6 + 6) = 0.0d0 + Urtot = Urtot + Ur + UadUrdf(Nppt3 + 1) = Ua + UadUrdf(Nppt3 + 2) = -alrij*Ur + UadUrdf(Nppt3 + 3) = -alaij*Ua + ELSE + PL1 = pi/(R2ij - R1ij) + PL2 = PL1*(Rij - R1ij) + fij = 0.5d0 + 0.5d0*dcos(PL2) + dfij = -0.5d0*PL1*dsin(PL2) + XYZRrefdf(Nppt6 + 5) = fij + XYZRrefdf(Nppt6 + 6) = dfij + Urtot = Urtot + fij*Ur + UadUrdf(Nppt3 + 1) = fij*Ua + UadUrdf(Nppt3 + 2) = (dfij - alrij*fij)*Ur + UadUrdf(Nppt3 + 3) = (dfij - alaij*fij)*Ua + END IF + END DO DO_J + END SUBROUTINE tersoff_subeniat_l + +! ************************************************************************************************** +!> \brief ... +!> \param i ... +!> \param Nmol ... +!> \param Npmax ... +!> \param NNmax ... +!> \param Kinds ... +!> \param Pn ... +!> \param Co_bcd ... +!> \param bcsq ... +!> \param dsq ... +!> \param h ... +!> \param XYZRrefdf ... +!> \param UadUrdf ... +!> \param F ... +!> \param Uatot ... +!> \param dkEij ... +!> \param lsta ... +!> \param lstb ... +!> \param nnbrx ... +! ************************************************************************************************** + SUBROUTINE tersoff_subfiat_l(i,Nmol,Npmax,NNmax,Kinds,Pn,Co_bcd,bcsq,dsq,h,XYZRrefdf,UadUrdf,F,Uatot,dkEij,lsta,lstb,nnbrx) + INTEGER :: i + INTEGER, INTENT(in) :: Nmol, Npmax, NNmax + INTEGER, DIMENSION(1:Nmol), INTENT(in) :: Kinds + REAL(8), DIMENSION(1:2), INTENT(in) :: Pn, Co_bcd, bcsq, dsq, h + REAL(8), DIMENSION(1:6*Npmax), INTENT(in) :: XYZRrefdf + REAL(8), DIMENSION(1:3*Npmax), INTENT(in) :: UadUrdf + REAL(8), DIMENSION(1:3*Nmol), INTENT(inout) :: F + REAL(8), INTENT(inout) :: Uatot + REAL(8), DIMENSION(1:3*NNmax) :: dkEij + INTEGER, INTENT(in) :: lsta(2, Nmol), nnbrx, lstb(nnbrx*Nmol) + + INTEGER :: ij, ijpt3, ijpt6, ik, ikpt6, Ipb, Ipe, & + Ipt3, Jpt3, Ki, Kpt3, Nkpt3 + REAL(8) :: bcsqi, Bij, Co1_dkEij, Co2_dkEij, Co_cdi, Co_dhcosi, Co_hcosi, Co_mb1, Co_mb2, & + Co_pa, COSijk, dfij, dfik, dFxi, dFxj, dFxk, dFyi, dFyj, dFyk, dFzi, dFzj, dFzk, dGi, & + djEij, dsqi, dXjEij2, dYjEij2, dZjEij2, Eij, fdG, fdGcos, fij, fik, Gi, hi, Pni, Rreij, & + Rreik, Ua, XRreij, XRreik, YRreij, YRreik, ZRreij, ZRreik + + Ipb = lsta(1, i) + Ipe = lsta(2, i) + + Ki = Kinds(i) + bcsqi = bcsq(Ki) + dsqi = dsq(Ki) + hi = h(Ki) + Pni = Pn(Ki) + + Co_cdi = Co_bcd(Ki) + + dFxi = 0.0d0 + dFyi = 0.0d0 + dFzi = 0.0d0 + + DO_J: DO ij = Ipb, Ipe, +1 + + IJpt3 = 3*(ij - 1) + IJpt6 = 6*(ij - 1) + + XRreij = XYZRrefdf(IJpt6 + 1) + YRreij = XYZRrefdf(IJpt6 + 2) + ZRreij = XYZRrefdf(IJpt6 + 3) + Rreij = XYZRrefdf(IJpt6 + 4) + fij = XYZRrefdf(IJpt6 + 5) + dfij = XYZRrefdf(IJpt6 + 6) + + Eij = 0.0d0 + djEij = 0.0d0 + dXjEij2 = 0.0d0 + dYjEij2 = 0.0d0 + dZjEij2 = 0.0d0 + + Nkpt3 = -3 + DO_K: DO ik = Ipb, Ipe, +1 + + Nkpt3 = Nkpt3 + 3 + + IKIJ: IF (ik /= ij) THEN + + IKpt6 = 6*(ik - 1) + + XRreik = XYZRrefdf(IKpt6 + 1) + YRreik = XYZRrefdf(IKpt6 + 2) + ZRreik = XYZRrefdf(IKpt6 + 3) + Rreik = XYZRrefdf(IKpt6 + 4) + fik = XYZRrefdf(IKpt6 + 5) + dfik = XYZRrefdf(IKpt6 + 6) + + COSijk = XRreij*XRreik + YRreij*YRreik + ZRreij*ZRreik + + Co_hcosi = hi - COSijk + Co_dhcosi = 1.0d0/(dsqi + Co_hcosi*Co_hcosi) + Gi = -bcsqi*Co_dhcosi + dGi = 2.0d0*Co_hcosi*Co_dhcosi*Gi + Gi = Gi + Co_cdi + + Eij = Eij + fik*Gi + + fdG = fik*dGi + fdGcos = fdG*COSijk + + djEij = djEij + fdGcos + + dXjEij2 = dXjEij2 + fdG*XRreik + dYjEij2 = dYjEij2 + fdG*YRreik + dZjEij2 = dZjEij2 + fdG*ZRreik + + Co1_dkEij = -dfik*Gi + fdGcos*Rreik + Co2_dkEij = -fdG*Rreik + + dkEij(Nkpt3 + 1) = Co1_dkEij*XRreik + Co2_dkEij*XRreij + dkEij(Nkpt3 + 2) = Co1_dkEij*YRreik + Co2_dkEij*YRreij + dkEij(Nkpt3 + 3) = Co1_dkEij*ZRreik + Co2_dkEij*ZRreij + + ELSE + dkEij(Nkpt3 + 1) = 0.0d0 + dkEij(Nkpt3 + 2) = 0.0d0 + dkEij(Nkpt3 + 3) = 0.0d0 + END IF IKIJ + + END DO DO_K + + Bij = 1.0d0 + Eij**Pni + Ua = UadUrdf(IJpt3 + 1)*Bij**(-0.5d0/Pni) + Uatot = Uatot + Ua + + Co_pa = UadUrdf(IJpt3 + 2) + UadUrdf(IJpt3 + 3)*Bij**(-0.5d0/Pni) + + CEij: IF (Nkpt3 > 0) THEN + + Co_mb1 = Ua*0.5d0*Eij**(Pni - 1.0d0)/Bij + Co_mb2 = Co_mb1*Rreij + + Nkpt3 = -3 + DO ik = Ipb, Ipe, +1 + + Nkpt3 = Nkpt3 + 3 + dFxk = Co_mb1*dkEij(Nkpt3 + 1) + dFyk = Co_mb1*dkEij(Nkpt3 + 2) + dFzk = Co_mb1*dkEij(Nkpt3 + 3) + + Kpt3 = 3*(lstb(ik) - 1) + F(Kpt3 + 1) = F(Kpt3 + 1) + dFxk + F(Kpt3 + 2) = F(Kpt3 + 2) + dFyk + F(Kpt3 + 3) = F(Kpt3 + 3) + dFzk + + dFxi = dFxi + dFxk + dFyi = dFyi + dFyk + dFzi = dFzi + dFzk + + END DO + + dFxj = Co_pa*XRreij + Co_mb2*(XRreij*djEij - dXjEij2) + dFyj = Co_pa*YRreij + Co_mb2*(YRreij*djEij - dYjEij2) + dFzj = Co_pa*ZRreij + Co_mb2*(ZRreij*djEij - dZjEij2) + + ELSE + + dFxj = Co_pa*XRreij + dFyj = Co_pa*YRreij + dFzj = Co_pa*ZRreij + + END IF CEij + + Jpt3 = 3*(lstb(ij) - 1) + + F(Jpt3 + 1) = F(Jpt3 + 1) + dFxj + F(Jpt3 + 2) = F(Jpt3 + 2) + dFyj + F(Jpt3 + 3) = F(Jpt3 + 3) + dFzj + + dFxi = dFxi + dFxj + dFyi = dFyi + dFyj + dFzi = dFzi + dFzj + + END DO DO_J + + Ipt3 = 3*(i - 1) + F(Ipt3 + 1) = F(Ipt3 + 1) - dFxi + F(Ipt3 + 2) = F(Ipt3 + 2) - dFyi + F(Ipt3 + 3) = F(Ipt3 + 3) - dFzi + END SUBROUTINE tersoff_subfiat_l + END MODULE eip_silicon diff --git a/src/force_env_methods.F b/src/force_env_methods.F index cc778c4043..01bbd43802 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -60,7 +60,9 @@ MODULE force_env_methods USE cp_units, ONLY: cp_unit_from_cp2k USE eip_environment_types, ONLY: eip_environment_type USE eip_silicon, ONLY: eip_bazant,& - eip_lenosky + eip_lenosky,& + eip_stillinger_weber,& + eip_tersoff USE embed_types, ONLY: embed_env_type,& opt_dmfet_pot_type,& opt_embed_pot_type @@ -87,7 +89,8 @@ MODULE force_env_methods USE grrm_utils, ONLY: write_grrm USE input_constants, ONLY: & debug_run, dfet, dmfet, mix_cdft, mix_coupled, mix_generic, mix_linear_combination, & - mix_minimum, mix_restrained, mixed_cdft_serial, use_bazant_eip, use_lenosky_eip + mix_minimum, mix_restrained, mixed_cdft_serial, use_bazant_eip, use_lenosky_eip, & + use_stillinger_weber_eip, use_tersoff_eip USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_retain,& section_vals_type,& @@ -271,11 +274,18 @@ CONTAINS e_gap = force_env%pwdft_env%energy%band_gap e_entropy = force_env%pwdft_env%energy%entropy CASE (use_eip_force) - IF (force_env%eip_env%eip_model == use_lenosky_eip) THEN + SELECT CASE (force_env%eip_env%eip_model) + CASE (use_lenosky_eip) CALL eip_lenosky(force_env%eip_env) - ELSE IF (force_env%eip_env%eip_model == use_bazant_eip) THEN + CASE (use_bazant_eip) CALL eip_bazant(force_env%eip_env) - END IF + CASE (use_stillinger_weber_eip) + CALL eip_stillinger_weber(force_env%eip_env) + CASE (use_tersoff_eip) + CALL eip_tersoff(force_env%eip_env) + CASE DEFAULT + CPABORT("Unknown EIP model.") + END SELECT CASE (use_qmmm) CALL qmmm_calc_energy_force(force_env%qmmm_env, & calculate_forces, energy_consistency, linres=linres_run) diff --git a/src/input_constants.F b/src/input_constants.F index b317be0fc1..0fd310ac8a 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -725,7 +725,9 @@ MODULE input_constants ! EIP models INTEGER, PARAMETER, PUBLIC :: use_bazant_eip = 1, & - use_lenosky_eip = 2 + use_lenosky_eip = 2, & + use_stillinger_weber_eip = 3, & + use_tersoff_eip = 4 ! ddapc restraint forms INTEGER, PARAMETER, PUBLIC :: do_ddapc_restraint = 773, & diff --git a/src/input_cp2k_eip.F b/src/input_cp2k_eip.F index d7f303712b..e1bddaa175 100644 --- a/src/input_cp2k_eip.F +++ b/src/input_cp2k_eip.F @@ -9,14 +9,22 @@ !> \brief Creates the EIP section of the input !> \par History !> 03.2006 created -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** MODULE input_cp2k_eip + USE bibliography, ONLY: Bazant1996,& + Bazant1997,& + Goedecker2002,& + Lenosky2000,& + Stillinger1985,& + Tersoff1988 USE cp_output_handling, ONLY: cp_print_key_section_create,& high_print_level,& medium_print_level USE input_constants, ONLY: use_bazant_eip,& - use_lenosky_eip + use_lenosky_eip,& + use_stillinger_weber_eip,& + use_tersoff_eip USE input_keyword_types, ONLY: keyword_create,& keyword_release,& keyword_type @@ -44,7 +52,7 @@ CONTAINS !> \param section the section to create !> \par History !> 03.2006 created -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE create_eip_section(section) TYPE(section_type), POINTER :: section @@ -58,19 +66,30 @@ CONTAINS CALL section_create(section, __LOCATION__, name="EIP", & description="This section contains all information to run an "// & "Empirical Interatomic Potential (EIP) calculation.", & - n_keywords=1, n_subsections=1, repeats=.FALSE.) + n_keywords=1, n_subsections=1, repeats=.FALSE., & + citations=[Bazant1996, Bazant1997, Goedecker2002, Lenosky2000, & + Stillinger1985, Tersoff1988]) NULLIFY (subsection, keyword) CALL keyword_create(keyword, __LOCATION__, name="EIP_MODEL", & - description="Selects the empirical interaction potential model", & + description="Selects the empirical interaction potential model. "// & + "EDIP is accepted as an alias of BAZANT and uses the identical "// & + "implementation. EIP is currently supported only for a single "// & + "MPI rank. BAZANT and LENOSKY retain OpenMP parallelization, "// & + "while STILLINGER_WEBER and TERSOFF currently run without "// & + "OpenMP parallelization.", & usage="EIP_MODEL BAZANT", type_of_var=enum_t, & n_var=1, repeats=.FALSE., variants=["EIP-MODEL"], & - enum_c_vals=s2a("BAZANT", "EDIP", "LENOSKY"), & - enum_i_vals=[use_bazant_eip, use_bazant_eip, use_lenosky_eip], & + enum_c_vals=s2a("BAZANT", "EDIP", "LENOSKY", & + "STILLINGER_WEBER", "TERSOFF"), & + enum_i_vals=[use_bazant_eip, use_bazant_eip, use_lenosky_eip, & + use_stillinger_weber_eip, use_tersoff_eip], & enum_desc=s2a("Bazant potentials", & "Environment-Dependent Interatomic Potential", & - "Lenosky potentials"), & + "Lenosky potentials", & + "Stillinger-Weber potentials", & + "Tersoff potentials"), & default_i_val=use_lenosky_eip) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -86,7 +105,7 @@ CONTAINS !> \param section the section to create !> \par History !> 03.2006 created -!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) +!> \author Thomas D. Kuehne (tkuehne@cp2k.org) ! ************************************************************************************************** SUBROUTINE create_eip_print_section(section) TYPE(section_type), POINTER :: section diff --git a/tests/Fist/regtest-1-1/Si_1000_bazant.inp b/tests/Fist/regtest-1-1/Si_1000_bazant.inp new file mode 100644 index 0000000000..1bb557181a --- /dev/null +++ b/tests/Fist/regtest-1-1/Si_1000_bazant.inp @@ -0,0 +1,30 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT Si_1000 + RUN_TYPE MD +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 10 + TEMPERATURE 300.0 + TIMESTEP 1.0 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD EIP + &EIP + EIP_MODEL BAZANT + &END EIP + &SUBSYS + &CELL + ABC 27.1550 27.1550 27.1550 + &END CELL + &TOPOLOGY + COORD_FILE_FORMAT XYZ + COORD_FILE_NAME ../sample_xyz/Si_1000.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/Fist/regtest-1-1/Si_1000_lenosky.inp b/tests/Fist/regtest-1-1/Si_1000_lenosky.inp new file mode 100644 index 0000000000..6bd9e4367f --- /dev/null +++ b/tests/Fist/regtest-1-1/Si_1000_lenosky.inp @@ -0,0 +1,30 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT Si_1000_Lenosky + RUN_TYPE MD +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 10 + TEMPERATURE 300.0 + TIMESTEP 1.0 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD EIP + &EIP + EIP_MODEL LENOSKY + &END EIP + &SUBSYS + &CELL + ABC 27.1550 27.1550 27.1550 + &END CELL + &TOPOLOGY + COORD_FILE_FORMAT XYZ + COORD_FILE_NAME ../sample_xyz/Si_1000.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/Fist/regtest-1-1/Si_1000_stillinger_weber.inp b/tests/Fist/regtest-1-1/Si_1000_stillinger_weber.inp new file mode 100644 index 0000000000..e72fb9890b --- /dev/null +++ b/tests/Fist/regtest-1-1/Si_1000_stillinger_weber.inp @@ -0,0 +1,30 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT Si_1000_SW + RUN_TYPE MD +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 10 + TEMPERATURE 300.0 + TIMESTEP 1.0 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD EIP + &EIP + EIP_MODEL STILLINGER_WEBER + &END EIP + &SUBSYS + &CELL + ABC 27.1550 27.1550 27.1550 + &END CELL + &TOPOLOGY + COORD_FILE_FORMAT XYZ + COORD_FILE_NAME ../sample_xyz/Si_1000.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/Fist/regtest-1-1/Si_1000_tersoff.inp b/tests/Fist/regtest-1-1/Si_1000_tersoff.inp new file mode 100644 index 0000000000..f03f37e7ec --- /dev/null +++ b/tests/Fist/regtest-1-1/Si_1000_tersoff.inp @@ -0,0 +1,30 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT Si_1000_Tersoff + RUN_TYPE MD +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 10 + TEMPERATURE 300.0 + TIMESTEP 1.0 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD EIP + &EIP + EIP_MODEL TERSOFF + &END EIP + &SUBSYS + &CELL + ABC 27.1550 27.1550 27.1550 + &END CELL + &TOPOLOGY + COORD_FILE_FORMAT XYZ + COORD_FILE_NAME ../sample_xyz/Si_1000.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/Fist/regtest-1-1/TEST_FILES.toml b/tests/Fist/regtest-1-1/TEST_FILES.toml index b673e72a1c..e38eaa9785 100644 --- a/tests/Fist/regtest-1-1/TEST_FILES.toml +++ b/tests/Fist/regtest-1-1/TEST_FILES.toml @@ -5,7 +5,10 @@ # for details see cp2k/tools/do_regtest # # EIP -"Si_1000.inp" = [{matcher="M002", tol=1.0E-14, ref=-0.170183162315E+03}] +"Si_1000_bazant.inp" = [{matcher="M002", tol=1.0E-14, ref=-0.170183162315E+03}] +"Si_1000_lenosky.inp" = [{matcher="M002", tol=1.0E-14, ref=-0.780798471651E+02}] +"Si_1000_stillinger_weber.inp" = [{matcher="M002", tol=1.0E-14, ref=-0.158634664291E+03}] +"Si_1000_tersoff.inp" = [{matcher="M002", tol=1.0E-14, ref=-0.169526339303E+03}] # CHARMM INTRA Forces Test "pot_input.inp" = [{matcher="M002", tol=1.0E-14, ref=0.427871243040E-01}] "pot_bond.inp" = [{matcher="M002", tol=1.0E-14, ref=0.183201063777E-01}] diff --git a/tests/Fist/regtest-1-1/Si_1000.inp b/tests/Fist/sample_xyz/Si_1000.xyz similarity index 99% rename from tests/Fist/regtest-1-1/Si_1000.inp rename to tests/Fist/sample_xyz/Si_1000.xyz index 2ab037439a..550fb89544 100644 --- a/tests/Fist/regtest-1-1/Si_1000.inp +++ b/tests/Fist/sample_xyz/Si_1000.xyz @@ -1,28 +1,5 @@ -&GLOBAL - PRINT_LEVEL LOW - PROJECT Si_1000 - RUN_TYPE MD -&END GLOBAL - -&MOTION - &MD - ENSEMBLE NVE - STEPS 10 - TEMPERATURE 300.0 - TIMESTEP 1.0 - &END MD -&END MOTION - -&FORCE_EVAL - METHOD EIP - &EIP - EIP_MODEL BAZANT - &END EIP - &SUBSYS - &CELL - ABC 27.1550 27.1550 27.1550 - &END CELL - &COORD +1000 +Si 1000 atom diamond supercell for EIP regtests Si 0.000000 0.000000 0.000000 Si 0.000000 2.715500 2.715500 Si 2.715500 0.000000 2.715500 @@ -1023,6 +1000,3 @@ Si 25.797250 25.797250 23.081750 Si 25.797250 23.081750 25.797250 Si 23.081750 25.797250 25.797250 - &END COORD - &END SUBSYS -&END FORCE_EVAL