diff --git a/lib/pace/CMakeLists.txt b/lib/pace/CMakeLists.txt
new file mode 100644
index 0000000000..2884c83cb9
--- /dev/null
+++ b/lib/pace/CMakeLists.txt
@@ -0,0 +1,19 @@
+cmake_minimum_required(VERSION 3.7) # CMake version check
+project(aceevaluator)
+set(CMAKE_CXX_STANDARD 11) # Enable c++11 standard
+
+
+set(PACE_EVALUATOR_PATH ${CMAKE_CURRENT_LIST_DIR}/src/USER-PACE)
+# message("CMakeLists.txt DEBUG: PACE_EVALUATOR_PATH=${PACE_EVALUATOR_PATH}")
+set(PACE_EVALUATOR_SRC_PATH ${PACE_EVALUATOR_PATH})
+
+FILE(GLOB PACE_EVALUATOR_SOURCE_FILES ${PACE_EVALUATOR_SRC_PATH}/*.cpp)
+list(FILTER PACE_EVALUATOR_SOURCE_FILES EXCLUDE REGEX ".*pair_pace.*")
+set(PACE_EVALUATOR_INCLUDE_DIR ${PACE_EVALUATOR_SRC_PATH})
+
+
+##### aceevaluator #####
+add_library(aceevaluator ${PACE_EVALUATOR_SOURCE_FILES})
+target_include_directories(aceevaluator PUBLIC ${PACE_EVALUATOR_INCLUDE_DIR})
+target_compile_options(aceevaluator PRIVATE -O3)
+set_target_properties(aceevaluator PROPERTIES OUTPUT_NAME pace${LAMMPS_MACHINE})
diff --git a/lib/pace/Install.py b/lib/pace/Install.py
index f87f5c14cb..640971e011 100644
--- a/lib/pace/Install.py
+++ b/lib/pace/Install.py
@@ -1 +1,111 @@
-# TODO
\ No newline at end of file
+# TODO#!/usr/bin/env python
+
+"""
+Install.py tool to download, compile, and setup the pace library
+used to automate the steps described in the README file in this dir
+"""
+
+from __future__ import print_function
+import sys, os, subprocess, shutil
+from argparse import ArgumentParser
+
+sys.path.append('..')
+from install_helpers import fullpath, geturl, checkmd5sum
+
+parser = ArgumentParser(prog='Install.py',
+ description="LAMMPS library build wrapper script")
+
+# settings
+
+thisdir = fullpath('.')
+version = "v.2021.2.3"
+
+# known checksums for different PACE versions. used to validate the download.
+checksums = { \
+ 'v.2021.2.3' : '9ebb087cba7e4ca041fde52f7e9e640c', \
+ }
+
+
+# help message
+
+HELP = """
+Syntax from src dir: make lib-pace args="-b"
+ or: make lib-pace args="-b -v version"
+Syntax from lib dir: python Install.py -b
+ or: python Install.py -b -v version
+
+Examples:
+
+make lib-pace args="-b" # install default version of PACE lib
+make lib-pace args="-b -v version" # install specified version of PACE lib
+
+
+"""
+
+pgroup = parser.add_mutually_exclusive_group()
+pgroup.add_argument("-b", "--build", action="store_true",
+ help="download and build base PACE library")
+parser.add_argument("-v", "--version", default=version, choices=checksums.keys(),
+ help="set version of PACE library to download and build (default: %s)" % version)
+parser.add_argument("-vv", "--verbose", action="store_true",
+ help="be more verbose about is happening while this script runs")
+
+args = parser.parse_args()
+
+# print help message and exit, if neither build nor path options are given
+if not args.build:
+ parser.print_help()
+ sys.exit(HELP)
+
+buildflag = args.build
+
+verboseflag = args.verbose
+version = args.version
+
+
+archive_extension = "tar.gz"
+url = "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/%s.%s" % (version, archive_extension)
+unarchived_folder_name = "lammps-user-pace-%s"%(version)
+
+# download PACE tarball, unpack, build PACE
+if buildflag:
+
+ # download entire tarball
+
+ print("Downloading pace tarball ...")
+ archive_filename = "%s.%s" % (version, archive_extension)
+ download_filename = "%s/%s" % (thisdir, archive_filename)
+ print("Downloading from ",url," to ",download_filename, end=" ")
+ geturl(url, download_filename)
+ print(" done")
+
+ # verify downloaded archive integrity via md5 checksum, if known.
+ if version in checksums:
+ if not checkmd5sum(checksums[version], archive_filename):
+ sys.exit("Checksum for pace library does not match")
+
+ print("Unpacking pace tarball ...")
+ src_folder = thisdir+"/src"
+ cmd = 'cd "%s"; rm -rf "%s"; tar -xvf %s; mv %s %s' % (thisdir, src_folder, archive_filename, unarchived_folder_name, src_folder)
+ subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
+
+ # configure
+
+ build_folder = "%s/build"%(thisdir)
+ print("Configuring libpace ...")
+ cmd = 'cd %s && mkdir build && cd build && cmake .. -DCMAKE_BUILD_TYPE=Release' % (thisdir)
+ txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+ if verboseflag: print(txt.decode("UTF-8"))
+
+ # build
+ print("Building libpace ...")
+ cmd = 'cd "%s" && make -j2 && cp libpace.a %s/' % (build_folder, thisdir)
+ txt = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
+ if verboseflag:
+ print(txt.decode("UTF-8"))
+
+# remove source files
+
+ print("Removing pace build files and archive ...")
+ cmd = 'rm %s; rm -rf %s' % (download_filename, build_folder)
+ subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
\ No newline at end of file
diff --git a/lib/pace/LICENSE b/lib/pace/LICENSE
deleted file mode 100644
index f288702d2f..0000000000
--- a/lib/pace/LICENSE
+++ /dev/null
@@ -1,674 +0,0 @@
- GNU GENERAL PUBLIC LICENSE
- Version 3, 29 June 2007
-
- Copyright (C) 2007 Free Software Foundation, Inc.
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
- Preamble
-
- The GNU General Public License is a free, copyleft license for
-software and other kinds of works.
-
- The licenses for most software and other practical works are designed
-to take away your freedom to share and change the works. By contrast,
-the GNU General Public License is intended to guarantee your freedom to
-share and change all versions of a program--to make sure it remains free
-software for all its users. We, the Free Software Foundation, use the
-GNU General Public License for most of our software; it applies also to
-any other work released this way by its authors. You can apply it to
-your programs, too.
-
- When we speak of free software, we are referring to freedom, not
-price. Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-them if you wish), that you receive source code or can get it if you
-want it, that you can change the software or use pieces of it in new
-free programs, and that you know you can do these things.
-
- To protect your rights, we need to prevent others from denying you
-these rights or asking you to surrender the rights. Therefore, you have
-certain responsibilities if you distribute copies of the software, or if
-you modify it: responsibilities to respect the freedom of others.
-
- For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must pass on to the recipients the same
-freedoms that you received. You must make sure that they, too, receive
-or can get the source code. And you must show them these terms so they
-know their rights.
-
- Developers that use the GNU GPL protect your rights with two steps:
-(1) assert copyright on the software, and (2) offer you this License
-giving you legal permission to copy, distribute and/or modify it.
-
- For the developers' and authors' protection, the GPL clearly explains
-that there is no warranty for this free software. For both users' and
-authors' sake, the GPL requires that modified versions be marked as
-changed, so that their problems will not be attributed erroneously to
-authors of previous versions.
-
- Some devices are designed to deny users access to install or run
-modified versions of the software inside them, although the manufacturer
-can do so. This is fundamentally incompatible with the aim of
-protecting users' freedom to change the software. The systematic
-pattern of such abuse occurs in the area of products for individuals to
-use, which is precisely where it is most unacceptable. Therefore, we
-have designed this version of the GPL to prohibit the practice for those
-products. If such problems arise substantially in other domains, we
-stand ready to extend this provision to those domains in future versions
-of the GPL, as needed to protect the freedom of users.
-
- Finally, every program is threatened constantly by software patents.
-States should not allow patents to restrict development and use of
-software on general-purpose computers, but in those that do, we wish to
-avoid the special danger that patents applied to a free program could
-make it effectively proprietary. To prevent this, the GPL assures that
-patents cannot be used to render the program non-free.
-
- The precise terms and conditions for copying, distribution and
-modification follow.
-
- TERMS AND CONDITIONS
-
- 0. Definitions.
-
- "This License" refers to version 3 of the GNU General Public License.
-
- "Copyright" also means copyright-like laws that apply to other kinds of
-works, such as semiconductor masks.
-
- "The Program" refers to any copyrightable work licensed under this
-License. Each licensee is addressed as "you". "Licensees" and
-"recipients" may be individuals or organizations.
-
- To "modify" a work means to copy from or adapt all or part of the work
-in a fashion requiring copyright permission, other than the making of an
-exact copy. The resulting work is called a "modified version" of the
-earlier work or a work "based on" the earlier work.
-
- A "covered work" means either the unmodified Program or a work based
-on the Program.
-
- To "propagate" a work means to do anything with it that, without
-permission, would make you directly or secondarily liable for
-infringement under applicable copyright law, except executing it on a
-computer or modifying a private copy. Propagation includes copying,
-distribution (with or without modification), making available to the
-public, and in some countries other activities as well.
-
- To "convey" a work means any kind of propagation that enables other
-parties to make or receive copies. Mere interaction with a user through
-a computer network, with no transfer of a copy, is not conveying.
-
- An interactive user interface displays "Appropriate Legal Notices"
-to the extent that it includes a convenient and prominently visible
-feature that (1) displays an appropriate copyright notice, and (2)
-tells the user that there is no warranty for the work (except to the
-extent that warranties are provided), that licensees may convey the
-work under this License, and how to view a copy of this License. If
-the interface presents a list of user commands or options, such as a
-menu, a prominent item in the list meets this criterion.
-
- 1. Source Code.
-
- The "source code" for a work means the preferred form of the work
-for making modifications to it. "Object code" means any non-source
-form of a work.
-
- A "Standard Interface" means an interface that either is an official
-standard defined by a recognized standards body, or, in the case of
-interfaces specified for a particular programming language, one that
-is widely used among developers working in that language.
-
- The "System Libraries" of an executable work include anything, other
-than the work as a whole, that (a) is included in the normal form of
-packaging a Major Component, but which is not part of that Major
-Component, and (b) serves only to enable use of the work with that
-Major Component, or to implement a Standard Interface for which an
-implementation is available to the public in source code form. A
-"Major Component", in this context, means a major essential component
-(kernel, window system, and so on) of the specific operating system
-(if any) on which the executable work runs, or a compiler used to
-produce the work, or an object code interpreter used to run it.
-
- The "Corresponding Source" for a work in object code form means all
-the source code needed to generate, install, and (for an executable
-work) run the object code and to modify the work, including scripts to
-control those activities. However, it does not include the work's
-System Libraries, or general-purpose tools or generally available free
-programs which are used unmodified in performing those activities but
-which are not part of the work. For example, Corresponding Source
-includes interface definition files associated with source files for
-the work, and the source code for shared libraries and dynamically
-linked subprograms that the work is specifically designed to require,
-such as by intimate data communication or control flow between those
-subprograms and other parts of the work.
-
- The Corresponding Source need not include anything that users
-can regenerate automatically from other parts of the Corresponding
-Source.
-
- The Corresponding Source for a work in source code form is that
-same work.
-
- 2. Basic Permissions.
-
- All rights granted under this License are granted for the term of
-copyright on the Program, and are irrevocable provided the stated
-conditions are met. This License explicitly affirms your unlimited
-permission to run the unmodified Program. The output from running a
-covered work is covered by this License only if the output, given its
-content, constitutes a covered work. This License acknowledges your
-rights of fair use or other equivalent, as provided by copyright law.
-
- You may make, run and propagate covered works that you do not
-convey, without conditions so long as your license otherwise remains
-in force. You may convey covered works to others for the sole purpose
-of having them make modifications exclusively for you, or provide you
-with facilities for running those works, provided that you comply with
-the terms of this License in conveying all material for which you do
-not control copyright. Those thus making or running the covered works
-for you must do so exclusively on your behalf, under your direction
-and control, on terms that prohibit them from making any copies of
-your copyrighted material outside their relationship with you.
-
- Conveying under any other circumstances is permitted solely under
-the conditions stated below. Sublicensing is not allowed; section 10
-makes it unnecessary.
-
- 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
-
- No covered work shall be deemed part of an effective technological
-measure under any applicable law fulfilling obligations under article
-11 of the WIPO copyright treaty adopted on 20 December 1996, or
-similar laws prohibiting or restricting circumvention of such
-measures.
-
- When you convey a covered work, you waive any legal power to forbid
-circumvention of technological measures to the extent such circumvention
-is effected by exercising rights under this License with respect to
-the covered work, and you disclaim any intention to limit operation or
-modification of the work as a means of enforcing, against the work's
-users, your or third parties' legal rights to forbid circumvention of
-technological measures.
-
- 4. Conveying Verbatim Copies.
-
- You may convey verbatim copies of the Program's source code as you
-receive it, in any medium, provided that you conspicuously and
-appropriately publish on each copy an appropriate copyright notice;
-keep intact all notices stating that this License and any
-non-permissive terms added in accord with section 7 apply to the code;
-keep intact all notices of the absence of any warranty; and give all
-recipients a copy of this License along with the Program.
-
- You may charge any price or no price for each copy that you convey,
-and you may offer support or warranty protection for a fee.
-
- 5. Conveying Modified Source Versions.
-
- You may convey a work based on the Program, or the modifications to
-produce it from the Program, in the form of source code under the
-terms of section 4, provided that you also meet all of these conditions:
-
- a) The work must carry prominent notices stating that you modified
- it, and giving a relevant date.
-
- b) The work must carry prominent notices stating that it is
- released under this License and any conditions added under section
- 7. This requirement modifies the requirement in section 4 to
- "keep intact all notices".
-
- c) You must license the entire work, as a whole, under this
- License to anyone who comes into possession of a copy. This
- License will therefore apply, along with any applicable section 7
- additional terms, to the whole of the work, and all its parts,
- regardless of how they are packaged. This License gives no
- permission to license the work in any other way, but it does not
- invalidate such permission if you have separately received it.
-
- d) If the work has interactive user interfaces, each must display
- Appropriate Legal Notices; however, if the Program has interactive
- interfaces that do not display Appropriate Legal Notices, your
- work need not make them do so.
-
- A compilation of a covered work with other separate and independent
-works, which are not by their nature extensions of the covered work,
-and which are not combined with it such as to form a larger program,
-in or on a volume of a storage or distribution medium, is called an
-"aggregate" if the compilation and its resulting copyright are not
-used to limit the access or legal rights of the compilation's users
-beyond what the individual works permit. Inclusion of a covered work
-in an aggregate does not cause this License to apply to the other
-parts of the aggregate.
-
- 6. Conveying Non-Source Forms.
-
- You may convey a covered work in object code form under the terms
-of sections 4 and 5, provided that you also convey the
-machine-readable Corresponding Source under the terms of this License,
-in one of these ways:
-
- a) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by the
- Corresponding Source fixed on a durable physical medium
- customarily used for software interchange.
-
- b) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by a
- written offer, valid for at least three years and valid for as
- long as you offer spare parts or customer support for that product
- model, to give anyone who possesses the object code either (1) a
- copy of the Corresponding Source for all the software in the
- product that is covered by this License, on a durable physical
- medium customarily used for software interchange, for a price no
- more than your reasonable cost of physically performing this
- conveying of source, or (2) access to copy the
- Corresponding Source from a network server at no charge.
-
- c) Convey individual copies of the object code with a copy of the
- written offer to provide the Corresponding Source. This
- alternative is allowed only occasionally and noncommercially, and
- only if you received the object code with such an offer, in accord
- with subsection 6b.
-
- d) Convey the object code by offering access from a designated
- place (gratis or for a charge), and offer equivalent access to the
- Corresponding Source in the same way through the same place at no
- further charge. You need not require recipients to copy the
- Corresponding Source along with the object code. If the place to
- copy the object code is a network server, the Corresponding Source
- may be on a different server (operated by you or a third party)
- that supports equivalent copying facilities, provided you maintain
- clear directions next to the object code saying where to find the
- Corresponding Source. Regardless of what server hosts the
- Corresponding Source, you remain obligated to ensure that it is
- available for as long as needed to satisfy these requirements.
-
- e) Convey the object code using peer-to-peer transmission, provided
- you inform other peers where the object code and Corresponding
- Source of the work are being offered to the general public at no
- charge under subsection 6d.
-
- A separable portion of the object code, whose source code is excluded
-from the Corresponding Source as a System Library, need not be
-included in conveying the object code work.
-
- A "User Product" is either (1) a "consumer product", which means any
-tangible personal property which is normally used for personal, family,
-or household purposes, or (2) anything designed or sold for incorporation
-into a dwelling. In determining whether a product is a consumer product,
-doubtful cases shall be resolved in favor of coverage. For a particular
-product received by a particular user, "normally used" refers to a
-typical or common use of that class of product, regardless of the status
-of the particular user or of the way in which the particular user
-actually uses, or expects or is expected to use, the product. A product
-is a consumer product regardless of whether the product has substantial
-commercial, industrial or non-consumer uses, unless such uses represent
-the only significant mode of use of the product.
-
- "Installation Information" for a User Product means any methods,
-procedures, authorization keys, or other information required to install
-and execute modified versions of a covered work in that User Product from
-a modified version of its Corresponding Source. The information must
-suffice to ensure that the continued functioning of the modified object
-code is in no case prevented or interfered with solely because
-modification has been made.
-
- If you convey an object code work under this section in, or with, or
-specifically for use in, a User Product, and the conveying occurs as
-part of a transaction in which the right of possession and use of the
-User Product is transferred to the recipient in perpetuity or for a
-fixed term (regardless of how the transaction is characterized), the
-Corresponding Source conveyed under this section must be accompanied
-by the Installation Information. But this requirement does not apply
-if neither you nor any third party retains the ability to install
-modified object code on the User Product (for example, the work has
-been installed in ROM).
-
- The requirement to provide Installation Information does not include a
-requirement to continue to provide support service, warranty, or updates
-for a work that has been modified or installed by the recipient, or for
-the User Product in which it has been modified or installed. Access to a
-network may be denied when the modification itself materially and
-adversely affects the operation of the network or violates the rules and
-protocols for communication across the network.
-
- Corresponding Source conveyed, and Installation Information provided,
-in accord with this section must be in a format that is publicly
-documented (and with an implementation available to the public in
-source code form), and must require no special password or key for
-unpacking, reading or copying.
-
- 7. Additional Terms.
-
- "Additional permissions" are terms that supplement the terms of this
-License by making exceptions from one or more of its conditions.
-Additional permissions that are applicable to the entire Program shall
-be treated as though they were included in this License, to the extent
-that they are valid under applicable law. If additional permissions
-apply only to part of the Program, that part may be used separately
-under those permissions, but the entire Program remains governed by
-this License without regard to the additional permissions.
-
- When you convey a copy of a covered work, you may at your option
-remove any additional permissions from that copy, or from any part of
-it. (Additional permissions may be written to require their own
-removal in certain cases when you modify the work.) You may place
-additional permissions on material, added by you to a covered work,
-for which you have or can give appropriate copyright permission.
-
- Notwithstanding any other provision of this License, for material you
-add to a covered work, you may (if authorized by the copyright holders of
-that material) supplement the terms of this License with terms:
-
- a) Disclaiming warranty or limiting liability differently from the
- terms of sections 15 and 16 of this License; or
-
- b) Requiring preservation of specified reasonable legal notices or
- author attributions in that material or in the Appropriate Legal
- Notices displayed by works containing it; or
-
- c) Prohibiting misrepresentation of the origin of that material, or
- requiring that modified versions of such material be marked in
- reasonable ways as different from the original version; or
-
- d) Limiting the use for publicity purposes of names of licensors or
- authors of the material; or
-
- e) Declining to grant rights under trademark law for use of some
- trade names, trademarks, or service marks; or
-
- f) Requiring indemnification of licensors and authors of that
- material by anyone who conveys the material (or modified versions of
- it) with contractual assumptions of liability to the recipient, for
- any liability that these contractual assumptions directly impose on
- those licensors and authors.
-
- All other non-permissive additional terms are considered "further
-restrictions" within the meaning of section 10. If the Program as you
-received it, or any part of it, contains a notice stating that it is
-governed by this License along with a term that is a further
-restriction, you may remove that term. If a license document contains
-a further restriction but permits relicensing or conveying under this
-License, you may add to a covered work material governed by the terms
-of that license document, provided that the further restriction does
-not survive such relicensing or conveying.
-
- If you add terms to a covered work in accord with this section, you
-must place, in the relevant source files, a statement of the
-additional terms that apply to those files, or a notice indicating
-where to find the applicable terms.
-
- Additional terms, permissive or non-permissive, may be stated in the
-form of a separately written license, or stated as exceptions;
-the above requirements apply either way.
-
- 8. Termination.
-
- You may not propagate or modify a covered work except as expressly
-provided under this License. Any attempt otherwise to propagate or
-modify it is void, and will automatically terminate your rights under
-this License (including any patent licenses granted under the third
-paragraph of section 11).
-
- However, if you cease all violation of this License, then your
-license from a particular copyright holder is reinstated (a)
-provisionally, unless and until the copyright holder explicitly and
-finally terminates your license, and (b) permanently, if the copyright
-holder fails to notify you of the violation by some reasonable means
-prior to 60 days after the cessation.
-
- Moreover, your license from a particular copyright holder is
-reinstated permanently if the copyright holder notifies you of the
-violation by some reasonable means, this is the first time you have
-received notice of violation of this License (for any work) from that
-copyright holder, and you cure the violation prior to 30 days after
-your receipt of the notice.
-
- Termination of your rights under this section does not terminate the
-licenses of parties who have received copies or rights from you under
-this License. If your rights have been terminated and not permanently
-reinstated, you do not qualify to receive new licenses for the same
-material under section 10.
-
- 9. Acceptance Not Required for Having Copies.
-
- You are not required to accept this License in order to receive or
-run a copy of the Program. Ancillary propagation of a covered work
-occurring solely as a consequence of using peer-to-peer transmission
-to receive a copy likewise does not require acceptance. However,
-nothing other than this License grants you permission to propagate or
-modify any covered work. These actions infringe copyright if you do
-not accept this License. Therefore, by modifying or propagating a
-covered work, you indicate your acceptance of this License to do so.
-
- 10. Automatic Licensing of Downstream Recipients.
-
- Each time you convey a covered work, the recipient automatically
-receives a license from the original licensors, to run, modify and
-propagate that work, subject to this License. You are not responsible
-for enforcing compliance by third parties with this License.
-
- An "entity transaction" is a transaction transferring control of an
-organization, or substantially all assets of one, or subdividing an
-organization, or merging organizations. If propagation of a covered
-work results from an entity transaction, each party to that
-transaction who receives a copy of the work also receives whatever
-licenses to the work the party's predecessor in interest had or could
-give under the previous paragraph, plus a right to possession of the
-Corresponding Source of the work from the predecessor in interest, if
-the predecessor has it or can get it with reasonable efforts.
-
- You may not impose any further restrictions on the exercise of the
-rights granted or affirmed under this License. For example, you may
-not impose a license fee, royalty, or other charge for exercise of
-rights granted under this License, and you may not initiate litigation
-(including a cross-claim or counterclaim in a lawsuit) alleging that
-any patent claim is infringed by making, using, selling, offering for
-sale, or importing the Program or any portion of it.
-
- 11. Patents.
-
- A "contributor" is a copyright holder who authorizes use under this
-License of the Program or a work on which the Program is based. The
-work thus licensed is called the contributor's "contributor version".
-
- A contributor's "essential patent claims" are all patent claims
-owned or controlled by the contributor, whether already acquired or
-hereafter acquired, that would be infringed by some manner, permitted
-by this License, of making, using, or selling its contributor version,
-but do not include claims that would be infringed only as a
-consequence of further modification of the contributor version. For
-purposes of this definition, "control" includes the right to grant
-patent sublicenses in a manner consistent with the requirements of
-this License.
-
- Each contributor grants you a non-exclusive, worldwide, royalty-free
-patent license under the contributor's essential patent claims, to
-make, use, sell, offer for sale, import and otherwise run, modify and
-propagate the contents of its contributor version.
-
- In the following three paragraphs, a "patent license" is any express
-agreement or commitment, however denominated, not to enforce a patent
-(such as an express permission to practice a patent or covenant not to
-sue for patent infringement). To "grant" such a patent license to a
-party means to make such an agreement or commitment not to enforce a
-patent against the party.
-
- If you convey a covered work, knowingly relying on a patent license,
-and the Corresponding Source of the work is not available for anyone
-to copy, free of charge and under the terms of this License, through a
-publicly available network server or other readily accessible means,
-then you must either (1) cause the Corresponding Source to be so
-available, or (2) arrange to deprive yourself of the benefit of the
-patent license for this particular work, or (3) arrange, in a manner
-consistent with the requirements of this License, to extend the patent
-license to downstream recipients. "Knowingly relying" means you have
-actual knowledge that, but for the patent license, your conveying the
-covered work in a country, or your recipient's use of the covered work
-in a country, would infringe one or more identifiable patents in that
-country that you have reason to believe are valid.
-
- If, pursuant to or in connection with a single transaction or
-arrangement, you convey, or propagate by procuring conveyance of, a
-covered work, and grant a patent license to some of the parties
-receiving the covered work authorizing them to use, propagate, modify
-or convey a specific copy of the covered work, then the patent license
-you grant is automatically extended to all recipients of the covered
-work and works based on it.
-
- A patent license is "discriminatory" if it does not include within
-the scope of its coverage, prohibits the exercise of, or is
-conditioned on the non-exercise of one or more of the rights that are
-specifically granted under this License. You may not convey a covered
-work if you are a party to an arrangement with a third party that is
-in the business of distributing software, under which you make payment
-to the third party based on the extent of your activity of conveying
-the work, and under which the third party grants, to any of the
-parties who would receive the covered work from you, a discriminatory
-patent license (a) in connection with copies of the covered work
-conveyed by you (or copies made from those copies), or (b) primarily
-for and in connection with specific products or compilations that
-contain the covered work, unless you entered into that arrangement,
-or that patent license was granted, prior to 28 March 2007.
-
- Nothing in this License shall be construed as excluding or limiting
-any implied license or other defenses to infringement that may
-otherwise be available to you under applicable patent law.
-
- 12. No Surrender of Others' Freedom.
-
- If conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License. If you cannot convey a
-covered work so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you may
-not convey it at all. For example, if you agree to terms that obligate you
-to collect a royalty for further conveying from those to whom you convey
-the Program, the only way you could satisfy both those terms and this
-License would be to refrain entirely from conveying the Program.
-
- 13. Use with the GNU Affero General Public License.
-
- Notwithstanding any other provision of this License, you have
-permission to link or combine any covered work with a work licensed
-under version 3 of the GNU Affero General Public License into a single
-combined work, and to convey the resulting work. The terms of this
-License will continue to apply to the part which is the covered work,
-but the special requirements of the GNU Affero General Public License,
-section 13, concerning interaction through a network will apply to the
-combination as such.
-
- 14. Revised Versions of this License.
-
- The Free Software Foundation may publish revised and/or new versions of
-the GNU General Public License from time to time. Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
- Each version is given a distinguishing version number. If the
-Program specifies that a certain numbered version of the GNU General
-Public License "or any later version" applies to it, you have the
-option of following the terms and conditions either of that numbered
-version or of any later version published by the Free Software
-Foundation. If the Program does not specify a version number of the
-GNU General Public License, you may choose any version ever published
-by the Free Software Foundation.
-
- If the Program specifies that a proxy can decide which future
-versions of the GNU General Public License can be used, that proxy's
-public statement of acceptance of a version permanently authorizes you
-to choose that version for the Program.
-
- Later license versions may give you additional or different
-permissions. However, no additional obligations are imposed on any
-author or copyright holder as a result of your choosing to follow a
-later version.
-
- 15. Disclaimer of Warranty.
-
- THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
-APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
-HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
-OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
-THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
-IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
-ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
-
- 16. Limitation of Liability.
-
- IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
-THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
-GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
-USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
-DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
-PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
-EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
-SUCH DAMAGES.
-
- 17. Interpretation of Sections 15 and 16.
-
- If the disclaimer of warranty and limitation of liability provided
-above cannot be given local legal effect according to their terms,
-reviewing courts shall apply local law that most closely approximates
-an absolute waiver of all civil liability in connection with the
-Program, unless a warranty or assumption of liability accompanies a
-copy of the Program in return for a fee.
-
- END OF TERMS AND CONDITIONS
-
- How to Apply These Terms to Your New Programs
-
- If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
- To do so, attach the following notices to the program. It is safest
-to attach them to the start of each source file to most effectively
-state the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
-
- Copyright (C)
-
- This program 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 3 of the License, or
- (at your option) any later version.
-
- This program 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 this program. If not, see .
-
-Also add information on how to contact you by electronic and paper mail.
-
- If the program does terminal interaction, make it output a short
-notice like this when it starts in an interactive mode:
-
- Copyright (C)
- This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
- This is free software, and you are welcome to redistribute it
- under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License. Of course, your program's commands
-might be different; for a GUI interface, you would use an "about box".
-
- You should also get your employer (if you work as a programmer) or school,
-if any, to sign a "copyright disclaimer" for the program, if necessary.
-For more information on this, and how to apply and follow the GNU GPL, see
-.
-
- The GNU General Public License does not permit incorporating your program
-into proprietary programs. If your program is a subroutine library, you
-may consider it more useful to permit linking proprietary applications with
-the library. If this is what you want to do, use the GNU Lesser General
-Public License instead of this License. But first, please read
-.
diff --git a/lib/pace/README b/lib/pace/README
index e69de29bb2..e9a1ab820a 100644
--- a/lib/pace/README
+++ b/lib/pace/README
@@ -0,0 +1,9 @@
+This directory contains files required to use the USER-PACE package.
+
+You can type "make lib-pace" from the src directory to see help on
+how to download and build this library via make commands, or you can
+do the same thing by typing "python Install.py" from within this
+directory.
+
+
+
diff --git a/lib/pace/ace_abstract_basis.cpp b/lib/pace/ace_abstract_basis.cpp
deleted file mode 100644
index 3d8afafdaf..0000000000
--- a/lib/pace/ace_abstract_basis.cpp
+++ /dev/null
@@ -1,184 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Lysogorskiy Yury on 28.04.2020.
-
-#include "ace_abstract_basis.h"
-
-////embedding function
-////case nemb = 1 only implementation
-////F = sign(x)*( ( 1 - exp(-(w*x)^3) )*abs(x)^m + ((1/w)^(m-1))*exp(-(w*x)^3)*abs(x) )
-//// !! no prefactor wpre
-void Fexp(DOUBLE_TYPE x, DOUBLE_TYPE m, DOUBLE_TYPE &F, DOUBLE_TYPE &DF) {
- DOUBLE_TYPE w = 1.e6;
- DOUBLE_TYPE eps = 1e-10;
-
- DOUBLE_TYPE lambda = pow(1.0 / w, m - 1.0);
- if (abs(x) > eps) {
- DOUBLE_TYPE g;
- DOUBLE_TYPE a = abs(x);
- DOUBLE_TYPE am = pow(a, m);
- DOUBLE_TYPE w3x3 = pow(w * a, 3);
- DOUBLE_TYPE sign_factor = (signbit(x) ? -1 : 1);
- if (w3x3 > 30.0)
- g = 0.0;
- else
- g = exp(-w3x3);
-
- DOUBLE_TYPE omg = 1.0 - g;
- F = sign_factor * (omg * am + lambda * g * a);
- DOUBLE_TYPE dg = -3.0 * w * w * w * a * a * g;
- DF = m * pow(a, m - 1.0) * omg - am * dg + lambda * dg * a + lambda * g;
- } else {
- F = lambda * x;
- DF = lambda;
- }
-}
-
-
-//Scaled-shifted embedding function
-//F = sign(x)*( ( 1 - exp(-(w*x)^3) )*abs(x)^m + ((1/w)^(m-1))*exp(-(w*x)^3)*abs(x) )
-// !! no prefactor wpre
-void FexpShiftedScaled(DOUBLE_TYPE rho, DOUBLE_TYPE mexp, DOUBLE_TYPE &F, DOUBLE_TYPE &DF) {
- DOUBLE_TYPE eps = 1e-10;
- DOUBLE_TYPE a, xoff, yoff, nx, exprho;
-
- if (abs(mexp - 1.0) < eps) {
- F = rho;
- DF = 1;
- } else {
- a = abs(rho);
- exprho = exp(-a);
- nx = 1. / mexp;
- xoff = pow(nx, (nx / (1.0 - nx))) * exprho;
- yoff = pow(nx, (1 / (1.0 - nx))) * exprho;
- DOUBLE_TYPE sign_factor = (signbit(rho) ? -1 : 1);
- F = sign_factor * (pow(xoff + a, mexp) - yoff);
- DF = yoff + mexp * (-xoff + 1.0) * pow(xoff + a, mexp - 1.);
- }
-}
-
-void ACEAbstractBasisSet::inner_cutoff(DOUBLE_TYPE rho_core, DOUBLE_TYPE rho_cut, DOUBLE_TYPE drho_cut,
- DOUBLE_TYPE &fcut, DOUBLE_TYPE &dfcut) {
-
- DOUBLE_TYPE rho_low = rho_cut - drho_cut;
- if (rho_core >= rho_cut) {
- fcut = 0;
- dfcut = 0;
- } else if (rho_core <= rho_low) {
- fcut = 1;
- dfcut = 0;
- } else {
- fcut = 0.5 * (1 + cos(M_PI * (rho_core - rho_low) / drho_cut));
- dfcut = -0.5 * sin(M_PI * (rho_core - rho_low) / drho_cut) * M_PI / drho_cut;
- }
-}
-
-void ACEAbstractBasisSet::FS_values_and_derivatives(Array1D &rhos, DOUBLE_TYPE &value,
- Array1D &derivatives, DENSITY_TYPE ndensity) {
- DOUBLE_TYPE F, DF = 0, wpre, mexp;
- for (int p = 0; p < ndensity; p++) {
- wpre = FS_parameters.at(p * ndensity + 0);
- mexp = FS_parameters.at(p * ndensity + 1);
- if (this->npoti == "FinnisSinclair")
- Fexp(rhos(p), mexp, F, DF);
- else if (this->npoti == "FinnisSinclairShiftedScaled")
- FexpShiftedScaled(rhos(p), mexp, F, DF);
- value += F * wpre; // * weight (wpre)
- derivatives(p) = DF * wpre;// * weight (wpre)
- }
-}
-
-void ACEAbstractBasisSet::_clean() {
-
- delete[] elements_name;
- elements_name = nullptr;
- delete radial_functions;
- radial_functions = nullptr;
-}
-
-ACEAbstractBasisSet::ACEAbstractBasisSet(const ACEAbstractBasisSet &other) {
- ACEAbstractBasisSet::_copy_scalar_memory(other);
- ACEAbstractBasisSet::_copy_dynamic_memory(other);
-}
-
-ACEAbstractBasisSet &ACEAbstractBasisSet::operator=(const ACEAbstractBasisSet &other) {
- if (this != &other) {
- // deallocate old memory
- ACEAbstractBasisSet::_clean();
- //copy scalar values
- ACEAbstractBasisSet::_copy_scalar_memory(other);
- //copy dynamic memory
- ACEAbstractBasisSet::_copy_dynamic_memory(other);
- }
- return *this;
-}
-
-ACEAbstractBasisSet::~ACEAbstractBasisSet() {
- ACEAbstractBasisSet::_clean();
-}
-
-void ACEAbstractBasisSet::_copy_scalar_memory(const ACEAbstractBasisSet &src) {
- deltaSplineBins = src.deltaSplineBins;
- FS_parameters = src.FS_parameters;
- npoti = src.npoti;
-
- nelements = src.nelements;
- rankmax = src.rankmax;
- ndensitymax = src.ndensitymax;
- nradbase = src.nradbase;
- lmax = src.lmax;
- nradmax = src.nradmax;
- cutoffmax = src.cutoffmax;
-
- spherical_harmonics = src.spherical_harmonics;
-
- rho_core_cutoffs = src.rho_core_cutoffs;
- drho_core_cutoffs = src.drho_core_cutoffs;
-
-
- E0vals = src.E0vals;
-}
-
-void ACEAbstractBasisSet::_copy_dynamic_memory(const ACEAbstractBasisSet &src) {//allocate new memory
- if (src.elements_name == nullptr)
- throw runtime_error("Could not copy ACEAbstractBasisSet::elements_name - array not initialized");
- elements_name = new string[nelements];
- //copy
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- elements_name[mu] = src.elements_name[mu];
- }
- radial_functions = src.radial_functions->clone();
-}
-
-SPECIES_TYPE ACEAbstractBasisSet::get_species_index_by_name(const string &elemname) {
- for (SPECIES_TYPE t = 0; t < nelements; t++) {
- if (this->elements_name[t] == elemname)
- return t;
- }
- return -1;
-}
\ No newline at end of file
diff --git a/lib/pace/ace_abstract_basis.h b/lib/pace/ace_abstract_basis.h
deleted file mode 100644
index 165ea9496f..0000000000
--- a/lib/pace/ace_abstract_basis.h
+++ /dev/null
@@ -1,169 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-
-// Created by Lysogorskiy Yury on 28.04.2020.
-
-#ifndef ACE_EVALUATOR_ACE_ABSTRACT_BASIS_H
-#define ACE_EVALUATOR_ACE_ABSTRACT_BASIS_H
-
-#include
-#include
-
-#include "ace_c_basisfunction.h"
-#include "ace_contigous_array.h"
-#include "ace_radial.h"
-#include "ace_spherical_cart.h"
-#include "ace_types.h"
-
-using namespace std;
-
-/**
- * Abstract basis set class
- */
-class ACEAbstractBasisSet {
-public:
- SPECIES_TYPE nelements = 0; ///< number of elements in basis set
- RANK_TYPE rankmax = 0; ///< maximum value of rank
- DENSITY_TYPE ndensitymax = 0; ///< maximum number of densities \f$ \rho^{(p)} \f$
- NS_TYPE nradbase = 0; ///< maximum number of radial \f$\textbf{basis}\f$ function \f$ g_{k}(r) \f$
- LS_TYPE lmax = 0; ///< \f$ l_\textrm{max} \f$ - maximum value of orbital moment \f$ l \f$
- NS_TYPE nradmax = 0; ///< maximum number \f$ n \f$ of radial function \f$ R_{nl}(r) \f$
- DOUBLE_TYPE cutoffmax = 0; ///< maximum value of cutoff distance among all species in basis set
- DOUBLE_TYPE deltaSplineBins = 0; ///< Spline interpolation density
-
- string npoti = "FinnisSinclair"; ///< FS and embedding function combination
-
- string *elements_name = nullptr; ///< Array of elements name for mapping from index (0..nelements-1) to element symbol (string)
-
- AbstractRadialBasis *radial_functions = nullptr; ///< object to work with radial functions
- ACECartesianSphericalHarmonics spherical_harmonics; ///< object to work with spherical harmonics in Cartesian representation
-
-
- Array1D rho_core_cutoffs; ///< energy-based inner cut-off
- Array1D drho_core_cutoffs; ///< decay of energy-based inner cut-off
-
- vector FS_parameters; ///< parameters for cluster functional, see Eq.(3) in implementation notes or Eq.(53) in PRB 99, 014104 (2019)
-
- // E0 values
- Array1D E0vals;
-
- /**
- * Default empty constructor
- */
- ACEAbstractBasisSet() = default;
-
- // copy constructor, operator= and destructor (see. Rule of Three)
-
- /**
- * Copy constructor (see. Rule of Three)
- * @param other
- */
- ACEAbstractBasisSet(const ACEAbstractBasisSet &other);
-
- /**
- * operator= (see. Rule of Three)
- * @param other
- * @return
- */
- ACEAbstractBasisSet &operator=(const ACEAbstractBasisSet &other);
-
- /**
- * virtual destructor (see. Rule of Three)
- */
- virtual ~ACEAbstractBasisSet();
-
- /**
- * Computing cluster functional \f$ F(\rho_i^{(1)}, \dots, \rho_i^{(P)}) \f$
- * and its derivatives \f$ (\partial F/\partial\rho_i^{(1)}, \dots, \partial F/\partial \rho_i^{(P)} ) \f$
- * @param rhos array with densities \f$ \rho^{(p)} \f$
- * @param value (out) return value of cluster functional
- * @param derivatives (out) array of derivatives \f$ (\partial F/\partial\rho_i^{(1)}, \dots, \partial F/\partial \rho_i^{(P)} ) \f$
- * @param ndensity number \f$ P \f$ of densities to use
- */
- void FS_values_and_derivatives(Array1D &rhos, DOUBLE_TYPE &value, Array1D &derivatives,
- DENSITY_TYPE ndensity);
-
- /**
- * Computing hard core pairwise repulsive potential \f$ f_{cut}(\rho_i^{(\textrm{core})})\f$ and its derivative,
- * see Eq.(29) of implementation notes
- * @param rho_core value of \f$ \rho_i^{(\textrm{core})} \f$
- * @param rho_cut \f$ \rho_{cut}^{\mu_i} \f$ value
- * @param drho_cut \f$ \Delta_{cut}^{\mu_i} \f$ value
- * @param fcut (out) return inner cutoff function
- * @param dfcut (out) return derivative of inner cutoff function
- */
- static void inner_cutoff(DOUBLE_TYPE rho_core, DOUBLE_TYPE rho_cut, DOUBLE_TYPE drho_cut, DOUBLE_TYPE &fcut,
- DOUBLE_TYPE &dfcut);
-
-
- /**
- * Virtual method to save potential to file
- * @param filename file name
- */
- virtual void save(const string &filename) = 0;
-
- /**
- * Virtual method to load potential from file
- * @param filename file name
- */
- virtual void load(const string filename) = 0;
-
- /**
- * Get the species index by its element name
- * @param elemname element name
- * @return species index
- */
- SPECIES_TYPE get_species_index_by_name(const string &elemname);
-
-
- // routines for copying and cleaning dynamic memory of the class (see. Rule of Three)
-
- /**
- * Routine for clean the dynamically allocated memory\n
- * IMPORTANT! It must be idempotent for safety.
- */
- virtual void _clean();
-
- /**
- * Copy dynamic memory from src. Must be override and extended in derived classes!
- * @param src source object to copy from
- */
- virtual void _copy_dynamic_memory(const ACEAbstractBasisSet &src);
-
- /**
- * Copy scalar values from src. Must be override and extended in derived classes!
- * @param src source object to copy from
- */
- virtual void _copy_scalar_memory(const ACEAbstractBasisSet &src);
-};
-
-void Fexp(DOUBLE_TYPE rho, DOUBLE_TYPE mexp, DOUBLE_TYPE &F, DOUBLE_TYPE &DF);
-
-void FexpShiftedScaled(DOUBLE_TYPE rho, DOUBLE_TYPE mexp, DOUBLE_TYPE &F, DOUBLE_TYPE &DF);
-
-#endif //ACE_EVALUATOR_ACE_ABSTRACT_BASIS_H
diff --git a/lib/pace/ace_array2dlm.h b/lib/pace/ace_array2dlm.h
deleted file mode 100644
index 2b38602bc7..0000000000
--- a/lib/pace/ace_array2dlm.h
+++ /dev/null
@@ -1,579 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-
-// Created by Yury Lysogorskiy on 11.01.20.
-
-
-#ifndef ACE_ARRAY2DLM_H
-#define ACE_ARRAY2DLM_H
-
-#include
-#include
-
-#include "ace_arraynd.h"
-#include "ace_contigous_array.h"
-#include "ace_types.h"
-
-using namespace std;
-
-/**
- * Contiguous array to organize values by \f$ (l,m) \f$ indiced (orbital moment and its projection).
- * Only \f$ l_\textrm{max}\f$ should be provided, \f$ m = -l, \dots,l \f$
- * for \f$ l = 0, \dots, l_\textrm{max}\f$
- * @tparam T type of values to store
- */
-template
-class Array2DLM : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
-
- LS_TYPE lmax = 0; ///< orbital dimension \f$ l_{max} \f$
-
- bool is_proxy = false; ///< flag to show, if object is owning the memory or just represent it (proxying)dimensions
-
-public:
- /**
- * Default empty constructor
- */
- Array2DLM() = default;
-
- /**
- * Parametrized constructor
- * @param lmax maximum value of \f$ l \f$
- * @param array_name name of the array
- */
- explicit Array2DLM(LS_TYPE lmax, string array_name = "Array2DLM") {
- init(lmax, array_name);
- }
-
-
- /**
- * Constructor to create slices-proxy array, i.e. to provide access to the memory, but not to own it.
- * @param lmax maximum value of \f$ l \f$
- * @param data_ptr pointer to original data
- * @param array_name name of the array
- */
- Array2DLM(LS_TYPE lmax, T *data_ptr, string array_name = "Array2DLM") {
- this->lmax = lmax;
- this->size = (lmax + 1) * (lmax + 1);
- this->data = data_ptr;
- this->array_name = array_name;
- is_proxy = true;
- };
-
- /**
- * Destructor
- */
- ~Array2DLM() {
- if (!is_proxy) {
- if (data != nullptr) delete[] data;
- }
- data = nullptr;
- }
-
- /**
- * Initialize array, allocate memory
- * @param lmax maximum value of l
- * @param array_name name of the array
- */
- void init(LS_TYPE lmax, string array_name = "Array2DLM") {
- if (is_proxy) {
- char s[1024];
- sprintf(s, "Could not re-initialize proxy-array %s\n", this->array_name.c_str());
- throw logic_error(s);
- }
- this->lmax = lmax;
- this->array_name = array_name;
- //for m = -l .. l
- if (size != (lmax + 1) * (lmax + 1)) {
- size = (lmax + 1) * (lmax + 1);
- if (data) delete[] data;
- data = new T[size]{};
- memset(data, 0.0, size * sizeof(T));
- } else {
- memset(data, 0, size * sizeof(T));
- }
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-/**
- * Check if indices (l,m) are within array
- */
- void check_indices(LS_TYPE l, MS_TYPE m) const {
-
- if ((l < 0) | (l > lmax)) {
- fprintf(stderr, "%s: Index l = %d out of range (0, %d)\n", array_name.c_str(), l, lmax);
- exit(EXIT_FAILURE);
- }
-
- if ((m < -l) | (m > l)) {
- fprintf(stderr, "%s: Index m = %d out of range (%d, %d)\n", array_name.c_str(), m, -l, l);
- exit(EXIT_FAILURE);
- }
- size_t ii = l * (l + 1) + m;
- if (ii >= size) {
- fprintf(stderr, "%s: index = %ld out of range %ld\n", array_name.c_str(), ii, size);
- exit(EXIT_FAILURE);
- }
- }
-#endif
-
- /**
- * Accessing the array value by index (l,m) for reading
- * @param l
- * @param m
- * @return array value
- */
- inline const T &operator()(LS_TYPE l, MS_TYPE m) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(l, m);
-#endif
- //l^2 + l + m
- return data[l * (l + 1) + m];
- }
-
- /**
- * Accessing the array value by index (l,m) for writing
- * @param l
- * @param m
- * @return array value
- */
- inline T &operator()(LS_TYPE l, MS_TYPE m) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(l, m);
-#endif
- //l^2 + l + m
- return data[l * (l + 1) + m];
- }
-
- /**
- * Convert array to STL vector> container
- * @return vector> container
- */
- vector> to_vector() const {
- vector> res;
- res.resize(lmax + 1);
-
- for (int i = 0; i < lmax + 1; i++) {
- res[i].resize(i + 1);
- for (int j = 0; j < i + 1; j++) {
- res[i][j] = operator()(i, j);
- }
- }
- return res;
- }
-};
-
-/**
- * Contiguous array to organize values by \f$ (i_0, l , m) \f$ indices.
- * Only \f$ d_{0}, l_\textrm{max}\f$ should be provided: \f$ m = -l, \dots,l \f$
- * for \f$ l = 0, \dots, l_\textrm{max}\f$
- * @tparam T type of values to store
- */
-template
-class Array3DLM : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
-
- LS_TYPE lmax = 0; ///< orbital dimension \f$ l_{max} \f$
-
-
- size_t dim[1] = {0}; ///< linear dimension \f$ d_{0} \f$
-
- size_t s[1] = {0}; ///< strides for linear dimensions
-
- Array1D *> _proxy_slices; ///< slices representation
-public:
- /**
- * Default empty constructor
- */
- Array3DLM() = default;
-
- /**
- * Parametrized constructor
- * @param array_name name of the array
- */
- Array3DLM(string array_name) {
- this->array_name = array_name;
- };
-
- /**
- * Parametrized constructor
- * @param d0 maximum value of \f$ i_0 \f$
- * @param lmax maximum value of \f$ l \f$
- * @param array_name name of the array
- */
- explicit Array3DLM(size_t d0, LS_TYPE lmax, string array_name = "Array3DLM") {
- init(d0, lmax, array_name);
- }
-
- /**
- * Initialize array and its slices
- * @param d0 maximum value of \f$ i_0 \f$
- * @param lmax maximum value of \f$ l \f$
- * @param array_name name of the array
- */
- void init(size_t d0, LS_TYPE lmax, string array_name = "Array3DLM") {
- this->array_name = array_name;
- this->lmax = lmax;
- dim[0] = d0;
- s[0] = lmax * lmax;
- if (size != s[0] * dim[0]) {
- size = s[0] * dim[0];
- if (data) delete[] data;
- data = new T[size]{};
- memset(data, 0, size * sizeof(T));
- } else {
- memset(data, 0, size * sizeof(T));
- }
-
- _proxy_slices.set_array_name(array_name + "_proxy");
- //arrange proxy-slices
- _clear_proxies();
- _proxy_slices.resize(dim[0]);
- for (size_t i0 = 0; i0 < dim[0]; ++i0) {
- _proxy_slices(i0) = new Array2DLM(this->lmax, &this->data[i0 * s[0]],
- array_name + "_slice");
- }
- }
-
- /**
- * Release pointers to slices
- */
- void _clear_proxies() {
- for (size_t i0 = 0; i0 < _proxy_slices.get_dim(0); ++i0) {
- delete _proxy_slices(i0);
- _proxy_slices(i0) = nullptr;
- }
- }
-
- /**
- * Destructor, clear proxies
- */
- ~Array3DLM() {
- _clear_proxies();
- }
-
- /**
- * Resize array to new dimensions
- * @param d0
- * @param lmax
- */
- void resize(size_t d0, LS_TYPE lmax) {
- _clear_proxies();
- init(d0, lmax, this->array_name);
- }
-
- /**
- * Get array dimensions
- * @param d dimension index
- * @return dimension along axis 'd'
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
- /**
- * Check if indices (i0, l,m) are within array
- */
- void check_indices(size_t i0, LS_TYPE l, MS_TYPE m) const {
- if ((l < 0) | (l > lmax)) {
- fprintf(stderr, "%s: Index l = %d out of range (0, %d)\n", array_name.c_str(), l, lmax);
- exit(EXIT_FAILURE);
- }
-
- if ((m < -l) | (m > l)) {
- fprintf(stderr, "%s: Index m = %d out of range (%d, %d)\n", array_name.c_str(), m, -l, l);
- exit(EXIT_FAILURE);
- }
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- fprintf(stderr, "%s: index i0 = %ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- exit(EXIT_FAILURE);
- }
-
- size_t ii = i0 * s[0] + l * (l + 1) + m;
- if (ii >= size) {
- fprintf(stderr, "%s: index = %ld out of range %ld\n", array_name.c_str(), ii, size);
- exit(EXIT_FAILURE);
- }
- }
-#endif
-
- /**
- * Accessing the array value by index (i0,l,m) for reading
- * @param i0
- * @param l
- * @param m
- * @return array value
- */
- inline const T &operator()(size_t i0, LS_TYPE l, MS_TYPE m) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, l, m);
-#endif
- return data[i0 * s[0] + l * (l + 1) + m];
- }
-
- /**
- * Accessing the array value by index (i0,l,m) for writing
- * @param i0
- * @param l
- * @param m
- * @return array value
- */
- inline T &operator()(size_t i0, LS_TYPE l, MS_TYPE m) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, l, m);
-#endif
- return data[i0 * s[0] + l * (l + 1) + m];
- }
-
- /**
- * Return proxy Array2DLM pointing to i0, l=0, m=0 to read
- * @param i0
- * @return proxy Array2DLM pointing to i0, l=0, m=0
- */
- inline const Array2DLM &operator()(size_t i0) const {
- return *_proxy_slices(i0);
- }
-
- /**
- * Return proxy Array2DLM pointing to i0, l=0, m=0 to write
- * @param i0
- * @return proxy Array2DLM pointing to i0, l=0, m=0
- */
- inline Array2DLM &operator()(size_t i0) {
- return *_proxy_slices(i0);
- }
-};
-
-
-/**
- * Contiguous array to organize values by \f$ (i_0, i_1, l , m) \f$ indices.
- * Only \f$ d_{0}, d_{1}, l_\textrm{max}\f$ should be provided: \f$ m = -l, \dots,l \f$
- * for \f$ l = 0, \dots, l_\textrm{max}\f$
- * @tparam T type of values to store
- */
-template
-class Array4DLM : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
-
- LS_TYPE lmax = 0; ///< orbital dimension \f$ l_{max} \f$
- size_t dim[2] = {0, 0}; ///< linear dimension \f$ d_{0}, d_{1} \f$
- size_t s[2] = {0, 0}; ///< strides for linear dimensions
-
- Array2D *> _proxy_slices; ///< slices representation
-public:
- /**
- * Default empty constructor
- */
- Array4DLM() = default;
-
- /**
- * Parametrized constructor
- * @param array_name name of the array
- */
- Array4DLM(string array_name) {
- this->array_name = array_name;
- };
-
- /**
- * Parametrized constructor
- * @param d0 maximum value of \f$ i_0 \f$
- * @param d1 maximum value of \f$ i_1 \f$
- * @param lmax maximum value of \f$ l \f$
- * @param array_name name of the array
- */
- explicit Array4DLM(size_t d0, size_t d1, LS_TYPE lmax, string array_name = "Array4DLM") {
- init(d0, d1, lmax, array_name);
- }
-
- /**
- * Initialize array, reallocate memory and its slices
- * @param d0 maximum value of \f$ i_0 \f$
- * @param d1 maximum value of \f$ i_1 \f$
- * @param lmax maximum value of \f$ l \f$
- * @param array_name name of the array
- */
- void init(size_t d0, size_t d1, LS_TYPE lmax, string array_name = "Array4DLM") {
- this->array_name = array_name;
- this->lmax = lmax;
- dim[1] = d1;
- dim[0] = d0;
- s[1] = lmax * lmax;
- s[0] = s[1] * dim[1];
- if (size != s[0] * dim[0]) {
- size = s[0] * dim[0];
- if (data) delete[] data;
- data = new T[size]{};
- memset(data, 0, size * sizeof(T));
- } else {
- memset(data, 0, size * sizeof(T));
- }
-
- _proxy_slices.set_array_name(array_name + "_proxy");
- //release old memory if there is any
- _clear_proxies();
- //arrange proxy-slices
- _proxy_slices.resize(dim[0], dim[1]);
- for (size_t i0 = 0; i0 < dim[0]; ++i0)
- for (size_t i1 = 0; i1 < dim[1]; ++i1) {
- _proxy_slices(i0, i1) = new Array2DLM(this->lmax, &this->data[i0 * s[0] + i1 * s[1]],
- array_name + "_slice");
- }
- }
-
- /**
- * Release pointers to slices
- */
- void _clear_proxies() {
-
- for (size_t i0 = 0; i0 < _proxy_slices.get_dim(0); ++i0)
- for (size_t i1 = 0; i1 < _proxy_slices.get_dim(1); ++i1) {
- delete _proxy_slices(i0, i1);
- _proxy_slices(i0, i1) = nullptr;
- }
- }
-
- /**
- * Destructor, clear proxies
- */
- ~Array4DLM() {
- _clear_proxies();
- }
-
- /**
- * Deallocate memory, reallocate with the new dimensions
- * @param d0
- * @param lmax
- */
- void resize(size_t d0, size_t d1, LS_TYPE lmax) {
- _clear_proxies();
- init(d0, d1, lmax, this->array_name);
- }
-
- /**
- * Get array dimensions
- * @param d dimension index
- * @return dimension along axis 'd'
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
- /**
- * Check if indices (i0, l,m) are within array
- */
- void check_indices(size_t i0, size_t i1, LS_TYPE l, MS_TYPE m) const {
- if ((l < 0) | (l > lmax)) {
- fprintf(stderr, "%s: Index l = %d out of range (0, %d)\n", array_name.c_str(), l, lmax);
- exit(EXIT_FAILURE);
- }
-
- if ((m < -l) | (m > l)) {
- fprintf(stderr, "%s: Index m = %d out of range (%d, %d)\n", array_name.c_str(), m, -l, l);
- exit(EXIT_FAILURE);
- }
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- fprintf(stderr, "%s: index i0 = %ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- exit(EXIT_FAILURE);
- }
-
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- fprintf(stderr, "%s: index i1 = %ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- exit(EXIT_FAILURE);
- }
-
- size_t ii = i0 * s[0] + i1 * s[1] + l * (l + 1) + m;
- if (ii >= size) {
- fprintf(stderr, "%s: index = %ld out of range %ld\n", array_name.c_str(), ii, size);
- exit(EXIT_FAILURE);
- }
- }
-#endif
-
- /**
- * Accessing the array value by index (i0,l,m) for reading
- * @param i0
- * @param i1
- * @param l
- * @param m
- * @return array value
- */
- inline const T &operator()(size_t i0, size_t i1, LS_TYPE l, MS_TYPE m) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, l, m);
-#endif
- return data[i0 * s[0] + i1 * s[1] + l * (l + 1) + m];
- }
-
- /**
- * Accessing the array value by index (i0,l,m) for writing
- * @param i0
- * @param i1
- * @param l
- * @param m
- * @return array value
- */
- inline T &operator()(size_t i0, size_t i1, LS_TYPE l, MS_TYPE m) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, l, m);
-#endif
- return data[i0 * s[0] + i1 * s[1] + l * (l + 1) + m];
- }
-
- /**
- * Return proxy Array2DLM pointing to i0, i1, l=0, m=0 to read
- * @param i0
- * @param i1
- * @return proxy Array2DLM pointing to i0, l=0, m=0
- */
- inline const Array2DLM &operator()(size_t i0, size_t i1) const {
- return *_proxy_slices(i0, i1);
- }
-
- /**
- * Return proxy Array2DLM pointing to i0, i1, l=0, m=0 to write
- * @param i0
- * @param i1
- * @return proxy Array2DLM pointing to i0, l=0, m=0
- */
- inline Array2DLM &operator()(size_t i0, size_t i1) {
- return *_proxy_slices(i0, i1);
- }
-};
-
-#endif //ACE_ARRAY2DLM_H
\ No newline at end of file
diff --git a/lib/pace/ace_arraynd.h b/lib/pace/ace_arraynd.h
deleted file mode 100644
index 1044b5654e..0000000000
--- a/lib/pace/ace_arraynd.h
+++ /dev/null
@@ -1,1949 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-
-//automatically generate source code
-
-#ifndef ACE_MULTIARRAY_H
-#define ACE_MULTIARRAY_H
-
-#include
-#include
-#include
-
-#include "ace_contigous_array.h"
-
-using namespace std;
-
-
-/**
- * Multidimensional (1 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array1D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[1] = {0}; ///< dimensions
- size_t s[1] = {0}; ///< strides
- int ndim = 1; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array1D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array1D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array1D(size_t d0, const string &array_name = "Array1D", T *new_data = nullptr) {
- init(d0, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0) {
-
- dim[0] = d0;
-
- s[0] = 1;
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, const string &array_name = "Array1D", T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0) {
- init(d0, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0) {
- //check data size consistency
- size_t new_size = d0;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0);
-#endif
-
- return data[i0];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0);
-#endif
-
- return data[i0];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array1D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector to_vector() const {
- vector res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0] = operator()(i0);
-
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector &vec) {
- size_t d0 = 0;
- d0 = vec.size();
-
-
- init(d0, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- operator()(i0) = vec.at(i0);
-
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array1D(const vector &vec, const string &array_name = "Array1D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array1D &operator=(const vector &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-/**
- * Multidimensional (2 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array2D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[2] = {0}; ///< dimensions
- size_t s[2] = {0}; ///< strides
- int ndim = 2; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array2D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array2D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array2D(size_t d0, size_t d1, const string &array_name = "Array2D", T *new_data = nullptr) {
- init(d0, d1, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0, size_t d1) {
-
- dim[0] = d0;
- dim[1] = d1;
-
- s[1] = 1;
- s[0] = s[1] * dim[1];
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, size_t d1, const string &array_name = "Array2D", T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0, d1);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0, size_t d1) {
- init(d0, d1, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0, size_t d1) {
- //check data size consistency
- size_t new_size = d0 * d1;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0, d1);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0, size_t i1) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- char buf[1024];
- sprintf(buf, "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0, size_t i1) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1);
-#endif
-
- return data[i0 * s[0] + i1];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0, size_t i1) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1);
-#endif
-
- return data[i0 * s[0] + i1];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array2D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector> to_vector() const {
- vector> res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0].resize(dim[1]);
-
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- res[i0][i1] = operator()(i0, i1);
-
- }
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector> &vec) {
- size_t d0 = 0;
- size_t d1 = 0;
- d0 = vec.size();
-
- if (d0 > 0) {
- d1 = vec.at(0).size();
-
-
- }
-
- init(d0, d1, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- if (vec.at(i0).size() != d1)
- throw std::invalid_argument("Vector size is not constant at dimension 1");
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- operator()(i0, i1) = vec.at(i0).at(i1);
-
- }
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array2D(const vector> &vec, const string &array_name = "Array2D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array2D &operator=(const vector> &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- /**
- * operator= to flatten vector container
- * @param vec container
- */
- Array2D &operator=(const vector &vec) {
- this->set_flatten_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-/**
- * Multidimensional (3 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array3D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[3] = {0}; ///< dimensions
- size_t s[3] = {0}; ///< strides
- int ndim = 3; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array3D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array3D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array3D(size_t d0, size_t d1, size_t d2, const string &array_name = "Array3D", T *new_data = nullptr) {
- init(d0, d1, d2, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0, size_t d1, size_t d2) {
-
- dim[0] = d0;
- dim[1] = d1;
- dim[2] = d2;
-
- s[2] = 1;
- s[1] = s[2] * dim[2];
- s[0] = s[1] * dim[1];
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, size_t d1, size_t d2, const string &array_name = "Array3D", T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0, d1, d2);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0, size_t d1, size_t d2) {
- init(d0, d1, d2, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0, size_t d1, size_t d2) {
- //check data size consistency
- size_t new_size = d0 * d1 * d2;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0, d1, d2);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0, size_t i1, size_t i2) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- char buf[1024];
- sprintf(buf, "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i2 < 0) | (i2 >= dim[2])) {
- char buf[1024];
- sprintf(buf, "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0, size_t i1, size_t i2) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0, size_t i1, size_t i2) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array3D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector>> to_vector() const {
- vector>> res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0].resize(dim[1]);
-
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- res[i0][i1].resize(dim[2]);
-
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- res[i0][i1][i2] = operator()(i0, i1, i2);
-
- }
- }
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector>> &vec) {
- size_t d0 = 0;
- size_t d1 = 0;
- size_t d2 = 0;
- d0 = vec.size();
-
- if (d0 > 0) {
- d1 = vec.at(0).size();
- if (d1 > 0) {
- d2 = vec.at(0).at(0).size();
-
-
- }
- }
-
- init(d0, d1, d2, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- if (vec.at(i0).size() != d1)
- throw std::invalid_argument("Vector size is not constant at dimension 1");
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- if (vec.at(i0).at(i1).size() != d2)
- throw std::invalid_argument("Vector size is not constant at dimension 2");
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- operator()(i0, i1, i2) = vec.at(i0).at(i1).at(i2);
-
- }
- }
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array3D(const vector>> &vec, const string &array_name = "Array3D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array3D &operator=(const vector>> &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- /**
- * operator= to flatten vector container
- * @param vec container
- */
- Array3D &operator=(const vector &vec) {
- this->set_flatten_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-/**
- * Multidimensional (4 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array4D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[4] = {0}; ///< dimensions
- size_t s[4] = {0}; ///< strides
- int ndim = 4; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array4D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array4D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array4D(size_t d0, size_t d1, size_t d2, size_t d3, const string &array_name = "Array4D", T *new_data = nullptr) {
- init(d0, d1, d2, d3, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0, size_t d1, size_t d2, size_t d3) {
-
- dim[0] = d0;
- dim[1] = d1;
- dim[2] = d2;
- dim[3] = d3;
-
- s[3] = 1;
- s[2] = s[3] * dim[3];
- s[1] = s[2] * dim[2];
- s[0] = s[1] * dim[1];
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, size_t d1, size_t d2, size_t d3, const string &array_name = "Array4D", T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0, d1, d2, d3);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0, size_t d1, size_t d2, size_t d3) {
- init(d0, d1, d2, d3, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0, size_t d1, size_t d2, size_t d3) {
- //check data size consistency
- size_t new_size = d0 * d1 * d2 * d3;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0, d1, d2, d3);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0, size_t i1, size_t i2, size_t i3) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- char buf[1024];
- sprintf(buf, "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i2 < 0) | (i2 >= dim[2])) {
- char buf[1024];
- sprintf(buf, "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i3 < 0) | (i3 >= dim[3])) {
- char buf[1024];
- sprintf(buf, "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0, size_t i1, size_t i2, size_t i3) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0, size_t i1, size_t i2, size_t i3) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array4D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector>>> to_vector() const {
- vector>>> res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0].resize(dim[1]);
-
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- res[i0][i1].resize(dim[2]);
-
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- res[i0][i1][i2].resize(dim[3]);
-
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- res[i0][i1][i2][i3] = operator()(i0, i1, i2, i3);
-
- }
- }
- }
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector>>> &vec) {
- size_t d0 = 0;
- size_t d1 = 0;
- size_t d2 = 0;
- size_t d3 = 0;
- d0 = vec.size();
-
- if (d0 > 0) {
- d1 = vec.at(0).size();
- if (d1 > 0) {
- d2 = vec.at(0).at(0).size();
- if (d2 > 0) {
- d3 = vec.at(0).at(0).at(0).size();
-
-
- }
- }
- }
-
- init(d0, d1, d2, d3, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- if (vec.at(i0).size() != d1)
- throw std::invalid_argument("Vector size is not constant at dimension 1");
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- if (vec.at(i0).at(i1).size() != d2)
- throw std::invalid_argument("Vector size is not constant at dimension 2");
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- if (vec.at(i0).at(i1).at(i2).size() != d3)
- throw std::invalid_argument("Vector size is not constant at dimension 3");
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- operator()(i0, i1, i2, i3) = vec.at(i0).at(i1).at(i2).at(i3);
-
- }
- }
- }
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array4D(const vector>>> &vec, const string &array_name = "Array4D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array4D &operator=(const vector>>> &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- /**
- * operator= to flatten vector container
- * @param vec container
- */
- Array4D &operator=(const vector &vec) {
- this->set_flatten_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-/**
- * Multidimensional (5 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array5D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[5] = {0}; ///< dimensions
- size_t s[5] = {0}; ///< strides
- int ndim = 5; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array5D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array5D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array5D(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, const string &array_name = "Array5D",
- T *new_data = nullptr) {
- init(d0, d1, d2, d3, d4, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4) {
-
- dim[0] = d0;
- dim[1] = d1;
- dim[2] = d2;
- dim[3] = d3;
- dim[4] = d4;
-
- s[4] = 1;
- s[3] = s[4] * dim[4];
- s[2] = s[3] * dim[3];
- s[1] = s[2] * dim[2];
- s[0] = s[1] * dim[1];
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, const string &array_name = "Array5D",
- T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0, d1, d2, d3, d4);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4) {
- init(d0, d1, d2, d3, d4, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4) {
- //check data size consistency
- size_t new_size = d0 * d1 * d2 * d3 * d4;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0, d1, d2, d3, d4);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- char buf[1024];
- sprintf(buf, "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i2 < 0) | (i2 >= dim[2])) {
- char buf[1024];
- sprintf(buf, "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i3 < 0) | (i3 >= dim[3])) {
- char buf[1024];
- sprintf(buf, "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i4 < 0) | (i4 >= dim[4])) {
- char buf[1024];
- sprintf(buf, "%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3, i4);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3 * s[3] + i4];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3, i4);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3 * s[3] + i4];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array5D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector>>>> to_vector() const {
- vector>>>> res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0].resize(dim[1]);
-
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- res[i0][i1].resize(dim[2]);
-
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- res[i0][i1][i2].resize(dim[3]);
-
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- res[i0][i1][i2][i3].resize(dim[4]);
-
-
- for (int i4 = 0; i4 < dim[4]; i4++) {
- res[i0][i1][i2][i3][i4] = operator()(i0, i1, i2, i3, i4);
-
- }
- }
- }
- }
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector>>>> &vec) {
- size_t d0 = 0;
- size_t d1 = 0;
- size_t d2 = 0;
- size_t d3 = 0;
- size_t d4 = 0;
- d0 = vec.size();
-
- if (d0 > 0) {
- d1 = vec.at(0).size();
- if (d1 > 0) {
- d2 = vec.at(0).at(0).size();
- if (d2 > 0) {
- d3 = vec.at(0).at(0).at(0).size();
- if (d3 > 0) {
- d4 = vec.at(0).at(0).at(0).at(0).size();
-
-
- }
- }
- }
- }
-
- init(d0, d1, d2, d3, d4, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- if (vec.at(i0).size() != d1)
- throw std::invalid_argument("Vector size is not constant at dimension 1");
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- if (vec.at(i0).at(i1).size() != d2)
- throw std::invalid_argument("Vector size is not constant at dimension 2");
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- if (vec.at(i0).at(i1).at(i2).size() != d3)
- throw std::invalid_argument("Vector size is not constant at dimension 3");
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- if (vec.at(i0).at(i1).at(i2).at(i3).size() != d4)
- throw std::invalid_argument("Vector size is not constant at dimension 4");
-
- for (int i4 = 0; i4 < dim[4]; i4++) {
- operator()(i0, i1, i2, i3, i4) = vec.at(i0).at(i1).at(i2).at(i3).at(i4);
-
- }
- }
- }
- }
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array5D(const vector>>>> &vec, const string &array_name = "Array5D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array5D &operator=(const vector>>>> &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- /**
- * operator= to flatten vector container
- * @param vec container
- */
- Array5D &operator=(const vector &vec) {
- this->set_flatten_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-/**
- * Multidimensional (6 - dimensional) array of type T with contiguous memory layout.
- * If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @tparam T data type
- */
-template
-class Array6D : public ContiguousArrayND {
- using ContiguousArrayND::array_name;
- using ContiguousArrayND::data;
- using ContiguousArrayND::size;
- using ContiguousArrayND::is_proxy_;
-
- size_t dim[6] = {0}; ///< dimensions
- size_t s[6] = {0}; ///< strides
- int ndim = 6; ///< number of dimensions
-public:
-
- /**
- * Default empty constructor
- */
- Array6D() = default;
-
- /**
- * Parametrized constructor with array name
- * @param array_name name of array (for error logging)
- */
- Array6D(const string &array_name) { this->array_name = array_name; }
-
- /**
- * Parametrized constructor
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- Array6D(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, size_t d5, const string &array_name = "Array6D",
- T *new_data = nullptr) {
- init(d0, d1, d2, d3, d4, d5, array_name, new_data);
- }
-
- /**
- * Setup the dimensions and strides of array
- * @param d0,... array sizes for different dimensions
- */
- void set_dimensions(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, size_t d5) {
-
- dim[0] = d0;
- dim[1] = d1;
- dim[2] = d2;
- dim[3] = d3;
- dim[4] = d4;
- dim[5] = d5;
-
- s[5] = 1;
- s[4] = s[5] * dim[5];
- s[3] = s[4] * dim[4];
- s[2] = s[3] * dim[3];
- s[1] = s[2] * dim[2];
- s[0] = s[1] * dim[1];
-
- size = s[0] * dim[0];
- };
-
- /**
- * Initialize array storage, dimensions and strides
- * @param d0,... array sizes for different dimensions
- * @param array_name string name of the array
- */
- void init(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, size_t d5, const string &array_name = "Array6D",
- T *new_data = nullptr) {
- this->array_name = array_name;
-
- size_t old_size = size;
- set_dimensions(d0, d1, d2, d3, d4, d5);
-
- bool new_is_proxy = (new_data != nullptr);
-
- if (new_is_proxy) {
- if (!is_proxy_) delete[] data;
- data = new_data;
- } else {
- if (size != old_size) {
- T *old_data = data; //preserve the pointer to the old data
- data = new T[size]; // allocate new data
- if (old_data != nullptr) { //
- size_t min_size = old_size < size ? old_size : size;
- memcpy(data, old_data, min_size * sizeof(T));
- if (!is_proxy_) delete[] old_data;
- }
- //memset(data, 0, size * sizeof(T));
- } else {
- //memset(data, 0, size * sizeof(T));
- }
- }
-
- is_proxy_ = new_is_proxy;
- }
-
- /**
- * Resize array
- * @param d0,... array sizes for different dimensions
- */
- void resize(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, size_t d5) {
- init(d0, d1, d2, d3, d4, d5, this->array_name);
- }
-
- /**
- * Reshape array without reset the data
- * @param d0,... array sizes for different dimensions
- */
- void reshape(size_t d0, size_t d1, size_t d2, size_t d3, size_t d4, size_t d5) {
- //check data size consistency
- size_t new_size = d0 * d1 * d2 * d3 * d4 * d5;
- if (new_size != size)
- throw invalid_argument("Couldn't reshape array when the size is not conserved");
- set_dimensions(d0, d1, d2, d3, d4, d5);
- }
-
- /**
- * Get array size in dimension "d"
- * @param d dimension index
- */
- size_t get_dim(int d) const {
- return dim[d];
- }
-
-#ifdef MULTIARRAY_INDICES_CHECK
-
- /**
- * Check indices validity. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- void check_indices(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4, size_t i5) const {
-
- if ((i0 < 0) | (i0 >= dim[0])) {
- char buf[1024];
- sprintf(buf, "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i1 < 0) | (i1 >= dim[1])) {
- char buf[1024];
- sprintf(buf, "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i2 < 0) | (i2 >= dim[2])) {
- char buf[1024];
- sprintf(buf, "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i3 < 0) | (i3 >= dim[3])) {
- char buf[1024];
- sprintf(buf, "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i4 < 0) | (i4 >= dim[4])) {
- char buf[1024];
- sprintf(buf, "%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4] - 1);
- throw std::out_of_range(buf);
- }
-
- if ((i5 < 0) | (i5 >= dim[5])) {
- char buf[1024];
- sprintf(buf, "%s: index i5=%ld out of range (0, %ld)\n", array_name.c_str(), i5, dim[5] - 1);
- throw std::out_of_range(buf);
- }
-
- }
-
-#endif
-
- /**
- * Array access operator() for reading. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline const T &operator()(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4, size_t i5) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3, i4, i5);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3 * s[3] + i4 * s[4] + i5];
-
- }
-
- /**
- * Array access operator() for writing. If preprocessor macro MULTIARRAY_INDICES_CHECK is defined, then the check of index will
- * be performed before accessing memory. By default this is turned off.
- * @param i0,... indices
- */
- inline T &operator()(size_t i0, size_t i1, size_t i2, size_t i3, size_t i4, size_t i5) {
-#ifdef MULTIARRAY_INDICES_CHECK
- check_indices(i0, i1, i2, i3, i4, i5);
-#endif
-
- return data[i0 * s[0] + i1 * s[1] + i2 * s[2] + i3 * s[3] + i4 * s[4] + i5];
-
- }
-
- /**
- * Array comparison operator
- * @param other
- */
- bool operator==(const Array6D &other) const {
- //compare dimensions
- for (int d = 0; d < ndim; d++) {
- if (this->dim[d] != other.dim[d])
- return false;
- }
- return ContiguousArrayND::operator==(other);
- }
-
- /**
- * Convert to nested vector> container
- * @return vector container
- */
- vector>>>>> to_vector() const {
- vector>>>>> res;
-
- res.resize(dim[0]);
-
- for (int i0 = 0; i0 < dim[0]; i0++) {
- res[i0].resize(dim[1]);
-
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- res[i0][i1].resize(dim[2]);
-
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- res[i0][i1][i2].resize(dim[3]);
-
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- res[i0][i1][i2][i3].resize(dim[4]);
-
-
- for (int i4 = 0; i4 < dim[4]; i4++) {
- res[i0][i1][i2][i3][i4].resize(dim[5]);
-
-
- for (int i5 = 0; i5 < dim[5]; i5++) {
- res[i0][i1][i2][i3][i4][i5] = operator()(i0, i1, i2, i3, i4, i5);
-
- }
- }
- }
- }
- }
- }
-
-
- return res;
- } // end to_vector()
-
-
- /**
- * Set values to vector> container
- * @param vec container
- */
- void set_vector(const vector>>>>> &vec) {
- size_t d0 = 0;
- size_t d1 = 0;
- size_t d2 = 0;
- size_t d3 = 0;
- size_t d4 = 0;
- size_t d5 = 0;
- d0 = vec.size();
-
- if (d0 > 0) {
- d1 = vec.at(0).size();
- if (d1 > 0) {
- d2 = vec.at(0).at(0).size();
- if (d2 > 0) {
- d3 = vec.at(0).at(0).at(0).size();
- if (d3 > 0) {
- d4 = vec.at(0).at(0).at(0).at(0).size();
- if (d4 > 0) {
- d5 = vec.at(0).at(0).at(0).at(0).at(0).size();
-
-
- }
- }
- }
- }
- }
-
- init(d0, d1, d2, d3, d4, d5, array_name);
- for (int i0 = 0; i0 < dim[0]; i0++) {
- if (vec.at(i0).size() != d1)
- throw std::invalid_argument("Vector size is not constant at dimension 1");
-
- for (int i1 = 0; i1 < dim[1]; i1++) {
- if (vec.at(i0).at(i1).size() != d2)
- throw std::invalid_argument("Vector size is not constant at dimension 2");
-
- for (int i2 = 0; i2 < dim[2]; i2++) {
- if (vec.at(i0).at(i1).at(i2).size() != d3)
- throw std::invalid_argument("Vector size is not constant at dimension 3");
-
- for (int i3 = 0; i3 < dim[3]; i3++) {
- if (vec.at(i0).at(i1).at(i2).at(i3).size() != d4)
- throw std::invalid_argument("Vector size is not constant at dimension 4");
-
- for (int i4 = 0; i4 < dim[4]; i4++) {
- if (vec.at(i0).at(i1).at(i2).at(i3).at(i4).size() != d5)
- throw std::invalid_argument("Vector size is not constant at dimension 5");
-
- for (int i5 = 0; i5 < dim[5]; i5++) {
- operator()(i0, i1, i2, i3, i4, i5) = vec.at(i0).at(i1).at(i2).at(i3).at(i4).at(i5);
-
- }
- }
- }
- }
- }
- }
-
- }
-
- /**
- * Parametrized constructor from vector> container
- * @param vec container
- * @param array_name array name
- */
- Array6D(const vector>>>>> &vec, const string &array_name = "Array6D") {
- this->set_vector(vec);
- this->array_name = array_name;
- }
-
- /**
- * operator= to vector> container
- * @param vec container
- */
- Array6D &operator=(const vector>>>>> &vec) {
- this->set_vector(vec);
- return *this;
- }
-
-
- /**
- * operator= to flatten vector container
- * @param vec container
- */
- Array6D &operator=(const vector &vec) {
- this->set_flatten_vector(vec);
- return *this;
- }
-
-
- vector get_shape() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = dim[d];
- return sh;
- }
-
- vector get_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d];
- return sh;
- }
-
- vector get_memory_strides() const {
- vector sh(ndim);
- for (int d = 0; d < ndim; d++)
- sh[d] = s[d] * sizeof(T);
- return sh;
- }
-};
-
-
-#endif //ACE_MULTIARRAY_H
diff --git a/lib/pace/ace_c_basis.cpp b/lib/pace/ace_c_basis.cpp
deleted file mode 100644
index d1c55700b7..0000000000
--- a/lib/pace/ace_c_basis.cpp
+++ /dev/null
@@ -1,980 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Yury Lysogorskiy on 1.04.20.
-
-#include "ace_c_basis.h"
-#include "ships_radial.h"
-
-using namespace std;
-
-ACECTildeBasisSet::ACECTildeBasisSet(string filename) {
- load(filename);
-}
-
-ACECTildeBasisSet::ACECTildeBasisSet(const ACECTildeBasisSet &other) {
- ACECTildeBasisSet::_copy_scalar_memory(other);
- ACECTildeBasisSet::_copy_dynamic_memory(other);
- pack_flatten_basis();
-}
-
-
-ACECTildeBasisSet &ACECTildeBasisSet::operator=(const ACECTildeBasisSet &other) {
- if (this != &other) {
- _clean();
- _copy_scalar_memory(other);
- _copy_dynamic_memory(other);
- pack_flatten_basis();
- }
- return *this;
-}
-
-
-ACECTildeBasisSet::~ACECTildeBasisSet() {
- ACECTildeBasisSet::_clean();
-}
-
-void ACECTildeBasisSet::_clean() {
- // call parent method
- ACEFlattenBasisSet::_clean();
- _clean_contiguous_arrays();
- _clean_basis_arrays();
-}
-
-void ACECTildeBasisSet::_copy_scalar_memory(const ACECTildeBasisSet &src) {
- ACEFlattenBasisSet::_copy_scalar_memory(src);
- num_ctilde_max = src.num_ctilde_max;
-}
-
-void ACECTildeBasisSet::_copy_dynamic_memory(const ACECTildeBasisSet &src) {//allocate new memory
- ACEFlattenBasisSet::_copy_dynamic_memory(src);
-
- if (src.basis_rank1 == nullptr)
- throw runtime_error("Could not copy ACECTildeBasisSet::basis_rank1 - array not initialized");
- if (src.basis == nullptr) throw runtime_error("Could not copy ACECTildeBasisSet::basis - array not initialized");
-
-
- basis_rank1 = new ACECTildeBasisFunction *[src.nelements];
- basis = new ACECTildeBasisFunction *[src.nelements];
-
- //copy basis arrays
- for (SPECIES_TYPE mu = 0; mu < src.nelements; ++mu) {
- basis_rank1[mu] = new ACECTildeBasisFunction[src.total_basis_size_rank1[mu]];
-
- for (size_t i = 0; i < src.total_basis_size_rank1[mu]; i++) {
- basis_rank1[mu][i] = src.basis_rank1[mu][i];
- }
-
-
- basis[mu] = new ACECTildeBasisFunction[src.total_basis_size[mu]];
- for (size_t i = 0; i < src.total_basis_size[mu]; i++) {
- basis[mu][i] = src.basis[mu][i];
- }
- }
- //DON"T COPY CONTIGUOUS ARRAY, REBUILD THEM
-}
-
-void ACECTildeBasisSet::_clean_contiguous_arrays() {
- ACEFlattenBasisSet::_clean_contiguous_arrays();
-
- delete[] full_c_tildes_rank1;
- full_c_tildes_rank1 = nullptr;
-
- delete[] full_c_tildes;
- full_c_tildes = nullptr;
-}
-
-//re-pack the constituent dynamic arrays of all basis functions in contiguous arrays
-void ACECTildeBasisSet::pack_flatten_basis() {
- compute_array_sizes(basis_rank1, basis);
-
- //1. clean contiguous arrays
- _clean_contiguous_arrays();
- //2. allocate contiguous arrays
- delete[] full_ns_rank1;
- full_ns_rank1 = new NS_TYPE[rank_array_total_size_rank1];
- delete[] full_ls_rank1;
- full_ls_rank1 = new NS_TYPE[rank_array_total_size_rank1];
- delete[] full_mus_rank1;
- full_mus_rank1 = new SPECIES_TYPE[rank_array_total_size_rank1];
- delete[] full_ms_rank1;
- full_ms_rank1 = new MS_TYPE[rank_array_total_size_rank1];
-
- delete[] full_c_tildes_rank1;
- full_c_tildes_rank1 = new DOUBLE_TYPE[coeff_array_total_size_rank1];
-
-
- delete[] full_ns;
- full_ns = new NS_TYPE[rank_array_total_size];
- delete[] full_ls;
- full_ls = new LS_TYPE[rank_array_total_size];
- delete[] full_mus;
- full_mus = new SPECIES_TYPE[rank_array_total_size];
- delete[] full_ms;
- full_ms = new MS_TYPE[ms_array_total_size];
-
- delete[] full_c_tildes;
- full_c_tildes = new DOUBLE_TYPE[coeff_array_total_size];
-
-
- //3. copy the values from private C_tilde_B_basis_function arrays to new contigous space
- //4. clean private memory
- //5. reassign private array pointers
-
- //r = 0, rank = 1
- size_t rank_array_ind_rank1 = 0;
- size_t coeff_array_ind_rank1 = 0;
- size_t ms_array_ind_rank1 = 0;
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- for (int func_ind_r1 = 0; func_ind_r1 < total_basis_size_rank1[mu]; ++func_ind_r1) {
- ACECTildeBasisFunction &func = basis_rank1[mu][func_ind_r1];
-
- //copy values ns from c_tilde_basis_function private memory to contigous memory part
- full_ns_rank1[rank_array_ind_rank1] = func.ns[0];
-
- //copy values ls from c_tilde_basis_function private memory to contigous memory part
- full_ls_rank1[rank_array_ind_rank1] = func.ls[0];
-
- //copy values mus from c_tilde_basis_function private memory to contigous memory part
- full_mus_rank1[rank_array_ind_rank1] = func.mus[0];
-
- //copy values ctildes from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_c_tildes_rank1[coeff_array_ind_rank1], func.ctildes,
- func.ndensity * sizeof(DOUBLE_TYPE));
-
-
- //copy values mus from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_ms_rank1[ms_array_ind_rank1], func.ms_combs,
- func.num_ms_combs *
- func.rank * sizeof(MS_TYPE));
-
- //release memory of each ACECTildeBasisFunction if it is not proxy
- func._clean();
-
- func.mus = &full_mus_rank1[rank_array_ind_rank1];
- func.ns = &full_ns_rank1[rank_array_ind_rank1];
- func.ls = &full_ls_rank1[rank_array_ind_rank1];
- func.ms_combs = &full_ms_rank1[ms_array_ind_rank1];
- func.ctildes = &full_c_tildes_rank1[coeff_array_ind_rank1];
- func.is_proxy = true;
-
- rank_array_ind_rank1 += func.rank;
- ms_array_ind_rank1 += func.rank *
- func.num_ms_combs;
- coeff_array_ind_rank1 += func.num_ms_combs * func.ndensity;
-
- }
- }
-
-
- //rank>1, r>0
- size_t rank_array_ind = 0;
- size_t coeff_array_ind = 0;
- size_t ms_array_ind = 0;
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- for (int func_ind = 0; func_ind < total_basis_size[mu]; ++func_ind) {
- ACECTildeBasisFunction &func = basis[mu][func_ind];
-
- //copy values mus from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_mus[rank_array_ind], func.mus,
- func.rank * sizeof(SPECIES_TYPE));
-
- //copy values ns from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_ns[rank_array_ind], func.ns,
- func.rank * sizeof(NS_TYPE));
- //copy values ls from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_ls[rank_array_ind], func.ls,
- func.rank * sizeof(LS_TYPE));
- //copy values mus from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_ms[ms_array_ind], func.ms_combs,
- func.num_ms_combs *
- func.rank * sizeof(MS_TYPE));
-
- //copy values ctildes from c_tilde_basis_function private memory to contigous memory part
- memcpy(&full_c_tildes[coeff_array_ind], func.ctildes,
- func.num_ms_combs * func.ndensity * sizeof(DOUBLE_TYPE));
-
-
- //release memory of each ACECTildeBasisFunction if it is not proxy
- func._clean();
-
- func.ns = &full_ns[rank_array_ind];
- func.ls = &full_ls[rank_array_ind];
- func.mus = &full_mus[rank_array_ind];
- func.ctildes = &full_c_tildes[coeff_array_ind];
- func.ms_combs = &full_ms[ms_array_ind];
- func.is_proxy = true;
-
- rank_array_ind += func.rank;
- ms_array_ind += func.rank *
- func.num_ms_combs;
- coeff_array_ind += func.num_ms_combs * func.ndensity;
- }
- }
-}
-
-void fwrite_c_tilde_b_basis_func(FILE *fptr, ACECTildeBasisFunction &func) {
- RANK_TYPE r;
- fprintf(fptr, "ctilde_basis_func: ");
- fprintf(fptr, "rank=%d ndens=%d mu0=%d ", func.rank, func.ndensity, func.mu0);
-
- fprintf(fptr, "mu=(");
- for (r = 0; r < func.rank; ++r)
- fprintf(fptr, " %d ", func.mus[r]);
- fprintf(fptr, ")\n");
-
- fprintf(fptr, "n=(");
- for (r = 0; r < func.rank; ++r)
- fprintf(fptr, " %d ", func.ns[r]);
- fprintf(fptr, ")\n");
-
- fprintf(fptr, "l=(");
- for (r = 0; r < func.rank; ++r)
- fprintf(fptr, " %d ", func.ls[r]);
- fprintf(fptr, ")\n");
-
- fprintf(fptr, "num_ms=%d\n", func.num_ms_combs);
-
- for (int m_ind = 0; m_ind < func.num_ms_combs; m_ind++) {
- fprintf(fptr, "<");
- for (r = 0; r < func.rank; ++r)
- fprintf(fptr, " %d ", func.ms_combs[m_ind * func.rank + r]);
- fprintf(fptr, ">: ");
- for (DENSITY_TYPE p = 0; p < func.ndensity; p++)
- fprintf(fptr, " %.18f ", func.ctildes[m_ind * func.ndensity + p]);
- fprintf(fptr, "\n");
- }
-
-}
-
-void ACECTildeBasisSet::save(const string &filename) {
- FILE *fptr;
- fptr = fopen(filename.c_str(), "w");
-
- fprintf(fptr, "nelements=%d\n", nelements);
-
- //elements mapping
- fprintf(fptr, "elements:");
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu)
- fprintf(fptr, " %s", elements_name[mu].c_str());
- fprintf(fptr, "\n\n");
-
- fprintf(fptr, "lmax=%d\n\n", lmax);
-
- fprintf(fptr, "embedding-function: %s\n", npoti.c_str());
-
- fprintf(fptr, "%ld FS parameters: ", FS_parameters.size());
- for (int i = 0; i < FS_parameters.size(); ++i) {
- fprintf(fptr, " %f", FS_parameters.at(i));
- }
- fprintf(fptr, "\n");
-
- //hard-core energy cutoff repulsion
- fprintf(fptr, "core energy-cutoff parameters: ");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- fprintf(fptr, "%.18f %.18f\n", rho_core_cutoffs(mu_i), drho_core_cutoffs(mu_i));
-
- // save E0 values
- fprintf(fptr, "E0:");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- fprintf(fptr, " %.18f", E0vals(mu_i));
- fprintf(fptr, "\n");
-
-
- fprintf(fptr, "\n");
-
-
- fprintf(fptr, "radbasename=%s\n", radial_functions->radbasename.c_str());
- fprintf(fptr, "nradbase=%d\n", nradbase);
- fprintf(fptr, "nradmax=%d\n", nradmax);
-
-
- fprintf(fptr, "cutoffmax=%f\n", cutoffmax);
-
- fprintf(fptr, "deltaSplineBins=%f\n", deltaSplineBins);
-
- //hard-core repulsion
- fprintf(fptr, "core repulsion parameters: ");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j)
- fprintf(fptr, "%.18f %.18f\n", radial_functions->prehc(mu_i, mu_j), radial_functions->lambdahc(mu_j, mu_j));
-
-
-
-
-
- //TODO: radial functions
- //radparameter
- fprintf(fptr, "radparameter=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j)
- fprintf(fptr, " %.18f", radial_functions->lambda(mu_i, mu_j));
- fprintf(fptr, "\n");
-
- fprintf(fptr, "cutoff=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j)
- fprintf(fptr, " %.18f", radial_functions->cut(mu_i, mu_j));
- fprintf(fptr, "\n");
-
- fprintf(fptr, "dcut=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j)
- fprintf(fptr, " %.18f", radial_functions->dcut(mu_i, mu_j));
- fprintf(fptr, "\n");
-
- fprintf(fptr, "crad=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
- for (NS_TYPE k = 0; k < nradbase; k++) {
- for (NS_TYPE n = 0; n < nradmax; n++) {
- for (LS_TYPE l = 0; l <= lmax; l++) {
- fprintf(fptr, " %.18f", radial_functions->crad(mu_i, mu_j, n, l, k));
- }
- fprintf(fptr, "\n");
- }
- }
- }
-
- fprintf(fptr, "\n");
-
- fprintf(fptr, "rankmax=%d\n", rankmax);
- fprintf(fptr, "ndensitymax=%d\n", ndensitymax);
- fprintf(fptr, "\n");
-
- //num_c_tilde_max
- fprintf(fptr, "num_c_tilde_max=%d\n", num_ctilde_max);
- fprintf(fptr, "num_ms_combinations_max=%d\n", num_ms_combinations_max);
-
-
- //write total_basis_size and total_basis_size_rank1
- fprintf(fptr, "total_basis_size_rank1: ");
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- fprintf(fptr, "%d ", total_basis_size_rank1[mu]);
- }
- fprintf(fptr, "\n");
-
- for (SPECIES_TYPE mu = 0; mu < nelements; mu++)
- for (SHORT_INT_TYPE func_ind = 0; func_ind < total_basis_size_rank1[mu]; ++func_ind)
- fwrite_c_tilde_b_basis_func(fptr, basis_rank1[mu][func_ind]);
-
- fprintf(fptr, "total_basis_size: ");
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- fprintf(fptr, "%d ", total_basis_size[mu]);
- }
- fprintf(fptr, "\n");
-
- for (SPECIES_TYPE mu = 0; mu < nelements; mu++)
- for (SHORT_INT_TYPE func_ind = 0; func_ind < total_basis_size[mu]; ++func_ind)
- fwrite_c_tilde_b_basis_func(fptr, basis[mu][func_ind]);
-
-
- fclose(fptr);
-}
-
-void fread_c_tilde_b_basis_func(FILE *fptr, ACECTildeBasisFunction &func) {
- RANK_TYPE r;
- int res;
- char buf[3][128];
-
- res = fscanf(fptr, " ctilde_basis_func: ");
-
- res = fscanf(fptr, "rank=%s ndens=%s mu0=%s ", buf[0], buf[1], buf[2]);
- if (res != 3)
- throw invalid_argument("Could not read C-tilde basis function");
-
- func.rank = (RANK_TYPE) stol(buf[0]);
- func.ndensity = (DENSITY_TYPE) stol(buf[1]);
- func.mu0 = (SPECIES_TYPE) stol(buf[2]);
-
- func.mus = new SPECIES_TYPE[func.rank];
- func.ns = new NS_TYPE[func.rank];
- func.ls = new LS_TYPE[func.rank];
-
- res = fscanf(fptr, " mu=(");
- for (r = 0; r < func.rank; ++r) {
- res = fscanf(fptr, "%s", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
- func.mus[r] = (SPECIES_TYPE) stol(buf[0]);
- }
- res = fscanf(fptr, " )"); // ")"
-
- res = fscanf(fptr, " n=("); // "n="
- for (r = 0; r < func.rank; ++r) {
- res = fscanf(fptr, "%s", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
-
- func.ns[r] = (NS_TYPE) stol(buf[0]);
- }
- res = fscanf(fptr, " )");
-
- res = fscanf(fptr, " l=(");
- for (r = 0; r < func.rank; ++r) {
- res = fscanf(fptr, "%s", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
- func.ls[r] = (NS_TYPE) stol(buf[0]);
- }
- res = fscanf(fptr, " )");
-
- res = fscanf(fptr, " num_ms=%s\n", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
-
- func.num_ms_combs = (SHORT_INT_TYPE) stoi(buf[0]);
-
- func.ms_combs = new MS_TYPE[func.rank * func.num_ms_combs];
- func.ctildes = new DOUBLE_TYPE[func.ndensity * func.num_ms_combs];
-
- for (int m_ind = 0; m_ind < func.num_ms_combs; m_ind++) {
- res = fscanf(fptr, " <");
- for (r = 0; r < func.rank; ++r) {
- res = fscanf(fptr, "%s", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
- func.ms_combs[m_ind * func.rank + r] = stoi(buf[0]);
- }
- res = fscanf(fptr, " >:");
- for (DENSITY_TYPE p = 0; p < func.ndensity; p++) {
- res = fscanf(fptr, "%s", buf[0]);
- if (res != 1)
- throw invalid_argument("Could not read C-tilde basis function");
- func.ctildes[m_ind * func.ndensity + p] = stod(buf[0]);
- }
- }
-}
-
-string
-format_error_message(const char *buffer, const string &filename, const string &var_name, const string &expected) {
- string err_message = "File '" + filename + "': couldn't read '" + var_name + "'";
- if (buffer)
- err_message = err_message + ", read:'" + buffer + "'";
- if (!expected.empty())
- err_message = err_message + ". Expected format: '" + expected + "'";
- return err_message;
-}
-
-void throw_error(const string &filename, const string &var_name, const string expected = "") {
- throw invalid_argument(format_error_message(nullptr, filename, var_name, expected));
-}
-
-DOUBLE_TYPE stod_err(const char *buffer, const string &filename, const string &var_name, const string expected = "") {
- try {
- return stod(buffer);
- } catch (invalid_argument &exc) {
- throw invalid_argument(format_error_message(buffer, filename, var_name, expected).c_str());
- }
-}
-
-int stoi_err(const char *buffer, const string &filename, const string &var_name, const string expected = "") {
- try {
- return stoi(buffer);
- } catch (invalid_argument &exc) {
- throw invalid_argument(format_error_message(buffer, filename, var_name, expected).c_str());
- }
-}
-
-long int stol_err(const char *buffer, const string &filename, const string &var_name, const string expected = "") {
- try {
- return stol(buffer);
- } catch (invalid_argument &exc) {
- throw invalid_argument(format_error_message(buffer, filename, var_name, expected).c_str());
- }
-}
-
-void ACECTildeBasisSet::load(const string filename) {
- int res;
- char buffer[1024], buffer2[1024];
- string radbasename = "ChebExpCos";
-
- FILE *fptr;
- fptr = fopen(filename.c_str(), "r");
- if (fptr == NULL)
- throw invalid_argument("Could not open file " + filename);
-
- //read number of elements
- res = fscanf(fptr, " nelements=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "nelements", "nelements=[number]");
- nelements = stoi_err(buffer, filename, "nelements", "nelements=[number]");
-
- //elements mapping
- elements_name = new string[nelements];
- res = fscanf(fptr, " elements:");
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "elements", "elements: Ele1 Ele2 ...");
- elements_name[mu] = buffer;
- }
-
- // load angular basis - only need spherical harmonics parameter
- res = fscanf(fptr, " lmax=%s\n", buffer);
- if (res != 1)
- throw_error(filename, "lmax", "lmax=[number]");
- lmax = stoi_err(buffer, filename, "lmax", "lmax=[number]");
- spherical_harmonics.init(lmax);
-
-
- // reading "embedding-function:"
- bool is_embedding_function_specified = false;
- res = fscanf(fptr, " embedding-function: %s", buffer);
- if (res == 0) {
- //throw_error(filename, "E0", " E0: E0-species1 E0-species2 ...");
- this->npoti = "FinnisSinclair"; // default values
- //printf("Warning! No embedding-function is specified, embedding-function: FinnisSinclair would be assumed\n");
- is_embedding_function_specified = false;
- } else {
- this->npoti = buffer;
- is_embedding_function_specified = true;
- }
- int parameters_size;
- res = fscanf(fptr, "%s FS parameters:", buffer);
- if (res != 1)
- throw_error(filename, "FS parameters size", "[number] FS parameters: par1 par2 ...");
- parameters_size = stoi_err(buffer, filename, "FS parameters size", "[number] FS parameters");
- FS_parameters.resize(parameters_size);
- for (int i = 0; i < FS_parameters.size(); ++i) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "FS parameter", "[number] FS parameters: [par1] [par2] ...");
- FS_parameters[i] = stod_err(buffer, filename, "FS parameter", "[number] FS parameters: [par1] [par2] ...");
- }
-
- if (!is_embedding_function_specified) {
- // assuming non-linear potential, and embedding function type is important
- for (int j = 1; j < parameters_size; j += 2)
- if (FS_parameters[j] != 1.0) { //if not ensure linearity of embedding function parameters
- printf("ERROR! Your potential is non-linear\n");
- printf("Please specify 'embedding-function: FinnisSinclair' or 'embedding-function: FinnisSinclairShiftedScaled' before 'FS parameters size' line\n");
- throw_error(filename, "embedding-function", "FinnisSinclair or FinnisSinclairShiftedScaled");
- }
- printf("Notice! No embedding-function is specified, but potential has linear embedding, default embedding-function: FinnisSinclair would be assumed\n");
- }
-
- //hard-core energy cutoff repulsion
- res = fscanf(fptr, " core energy-cutoff parameters:");
- if (res != 0)
- throw_error(filename, "core energy-cutoff parameters", "core energy-cutoff parameters: [par1] [par2]");
-
- rho_core_cutoffs.init(nelements, "rho_core_cutoffs");
- drho_core_cutoffs.init(nelements, "drho_core_cutoffs");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i) {
- res = fscanf(fptr, "%s %s", buffer, buffer2);
- if (res != 2)
- throw_error(filename, "core energy-cutoff parameters",
- "core energy-cutoff parameters: [rho_core_cut] [drho_core_cutoff] ...");
- rho_core_cutoffs(mu_i) = stod_err(buffer, filename, "rho core cutoff",
- "core energy-cutoff parameters: [rho_core_cut] [drho_core_cutoff] ...");
- drho_core_cutoffs(mu_i) = stod_err(buffer2, filename, "drho_core_cutoff",
- "core energy-cutoff parameters: [rho_core_cut] [drho_core_cutoff] ...");
- }
-
- // atom energy shift E0 (energy of isolated atom)
- E0vals.init(nelements);
-
- // reading "E0:"
- res = fscanf(fptr, " E0: %s", buffer);
- if (res == 0) {
- //throw_error(filename, "E0", " E0: E0-species1 E0-species2 ...");
- E0vals.fill(0.0);
- } else {
- double E0 = atof(buffer);
- E0vals(0) = E0;
-
- for (SPECIES_TYPE mu_i = 1; mu_i < nelements; ++mu_i) {
- res = fscanf(fptr, " %lf", &E0);
- if (res != 1)
- throw_error(filename, "E0", " couldn't read one of the E0 values");
- E0vals(mu_i) = E0;
- }
- res = fscanf(fptr, "\n");
- if (res != 0)
- printf("file %s : format seems broken near E0; trying to continue...\n", filename.c_str());
- }
-
- // check which radial basis we need to load
- res = fscanf(fptr, " radbasename=%s\n", buffer);
- if (res != 1) {
- throw_error(filename, "radbasename", "rabbasename=ChebExpCos|ChebPow|ACE.jl.Basic");
- } else {
- radbasename = buffer;
- }
-
-// printf("radbasename = `%s`\n", radbasename.c_str());
- if (radbasename == "ChebExpCos" | radbasename == "ChebPow") {
- _load_radial_ACERadial(fptr, filename, radbasename);
- } else if (radbasename == "ACE.jl.Basic") {
- _load_radial_SHIPsBasic(fptr, filename, radbasename);
- } else {
- throw invalid_argument(
- ("File '" + filename + "': I don't know how to read radbasename = " + radbasename).c_str());
- }
-
- res = fscanf(fptr, " rankmax=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "rankmax", "rankmax=[number]");
- rankmax = stoi_err(buffer, filename, "rankmax", "rankmax=[number]");
-
- res = fscanf(fptr, " ndensitymax=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "ndensitymax", "ndensitymax=[number]");
- ndensitymax = stoi_err(buffer, filename, "ndensitymax", "ndensitymax=[number]");
-
- // read the list of correlations to be put into the basis
- //num_c_tilde_max
- res = fscanf(fptr, " num_c_tilde_max=");
- res = fscanf(fptr, "%s\n", buffer);
- if (res != 1)
- throw_error(filename, "num_c_tilde_max", "num_c_tilde_max=[number]");
- num_ctilde_max = stol_err(buffer, filename, "num_c_tilde_max", "num_c_tilde_max=[number]");
-
- res = fscanf(fptr, " num_ms_combinations_max=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "num_ms_combinations_max", "num_ms_combinations_max=[number]");
-// throw invalid_argument(("File '" + filename + "': couldn't read num_ms_combinations_max").c_str());
- num_ms_combinations_max = stol_err(buffer, filename, "num_ms_combinations_max", "num_ms_combinations_max=[number]");
-
- //read total_basis_size_rank1
- total_basis_size_rank1 = new SHORT_INT_TYPE[nelements];
- basis_rank1 = new ACECTildeBasisFunction *[nelements];
- res = fscanf(fptr, " total_basis_size_rank1: ");
-
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "total_basis_size_rank1", "total_basis_size_rank1: [size_ele1] [size_ele2] ...");
-// throw invalid_argument(("File '" + filename + "': couldn't read total_basis_size_rank1").c_str());
- total_basis_size_rank1[mu] = stoi_err(buffer, filename, "total_basis_size_rank1",
- "total_basis_size_rank1: [size_ele1] [size_ele2] ...");
- basis_rank1[mu] = new ACECTildeBasisFunction[total_basis_size_rank1[mu]];
- }
- for (SPECIES_TYPE mu = 0; mu < nelements; mu++)
- for (SHORT_INT_TYPE func_ind = 0; func_ind < total_basis_size_rank1[mu]; ++func_ind) {
- fread_c_tilde_b_basis_func(fptr, basis_rank1[mu][func_ind]);
- }
-
- //read total_basis_size
- res = fscanf(fptr, " total_basis_size: ");
- total_basis_size = new SHORT_INT_TYPE[nelements];
- basis = new ACECTildeBasisFunction *[nelements];
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "total_basis_size", "total_basis_size: [size_ele1] [size_ele2] ...");
- total_basis_size[mu] = stoi_err(buffer, filename, "total_basis_size",
- "total_basis_size: [size_ele1] [size_ele2] ...");
- basis[mu] = new ACECTildeBasisFunction[total_basis_size[mu]];
- }
- for (SPECIES_TYPE mu = 0; mu < nelements; mu++)
- for (SHORT_INT_TYPE func_ind = 0; func_ind < total_basis_size[mu]; ++func_ind) {
- fread_c_tilde_b_basis_func(fptr, basis[mu][func_ind]);
- }
-
- fclose(fptr);
-
-// radial_functions->radbasename = radbasename;
- radial_functions->setuplookupRadspline();
- pack_flatten_basis();
-}
-
-void ACECTildeBasisSet::compute_array_sizes(ACECTildeBasisFunction **basis_rank1, ACECTildeBasisFunction **basis) {
- //compute arrays sizes
- rank_array_total_size_rank1 = 0;
- //ms_array_total_size_rank1 = rank_array_total_size_rank1;
- coeff_array_total_size_rank1 = 0;
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- if (total_basis_size_rank1[mu] > 0) {
- rank_array_total_size_rank1 += total_basis_size_rank1[mu];
-
- ACEAbstractBasisFunction &func = basis_rank1[mu][0];//TODO: get total density instead of density from first function
- coeff_array_total_size_rank1 += total_basis_size_rank1[mu] * func.ndensity;
- }
- }
-
- rank_array_total_size = 0;
- coeff_array_total_size = 0;
-
- ms_array_total_size = 0;
- max_dB_array_size = 0;
-
-
- max_B_array_size = 0;
-
- size_t cur_ms_size = 0;
- size_t cur_ms_rank_size = 0;
-
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
-
- cur_ms_size = 0;
- cur_ms_rank_size = 0;
- for (int func_ind = 0; func_ind < total_basis_size[mu]; ++func_ind) {
- auto &func = basis[mu][func_ind];
- rank_array_total_size += func.rank;
- ms_array_total_size += func.rank * func.num_ms_combs;
- coeff_array_total_size += func.ndensity * func.num_ms_combs;
-
- cur_ms_size += func.num_ms_combs;
- cur_ms_rank_size += func.rank * func.num_ms_combs;
- }
-
- if (cur_ms_size > max_B_array_size)
- max_B_array_size = cur_ms_size;
-
- if (cur_ms_rank_size > max_dB_array_size)
- max_dB_array_size = cur_ms_rank_size;
- }
-}
-
-void ACECTildeBasisSet::_clean_basis_arrays() {
- if (basis_rank1 != nullptr)
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- delete[] basis_rank1[mu];
- basis_rank1[mu] = nullptr;
- }
-
- if (basis != nullptr)
- for (SPECIES_TYPE mu = 0; mu < nelements; ++mu) {
- delete[] basis[mu];
- basis[mu] = nullptr;
- }
- delete[] basis;
- basis = nullptr;
-
- delete[] basis_rank1;
- basis_rank1 = nullptr;
-}
-
-//pack into 1D array with all basis functions
-void ACECTildeBasisSet::flatten_basis(C_tilde_full_basis_vector2d &mu0_ctilde_basis_vector) {
-
- _clean_basis_arrays();
- basis_rank1 = new ACECTildeBasisFunction *[nelements];
- basis = new ACECTildeBasisFunction *[nelements];
-
- delete[] total_basis_size_rank1;
- delete[] total_basis_size;
- total_basis_size_rank1 = new SHORT_INT_TYPE[nelements];
- total_basis_size = new SHORT_INT_TYPE[nelements];
-
-
- size_t tot_size_rank1 = 0;
- size_t tot_size = 0;
-
- for (SPECIES_TYPE mu = 0; mu < this->nelements; ++mu) {
- tot_size = 0;
- tot_size_rank1 = 0;
-
- for (auto &func: mu0_ctilde_basis_vector[mu]) {
- if (func.rank == 1) tot_size_rank1 += 1;
- else tot_size += 1;
- }
-
- total_basis_size_rank1[mu] = tot_size_rank1;
- basis_rank1[mu] = new ACECTildeBasisFunction[tot_size_rank1];
-
- total_basis_size[mu] = tot_size;
- basis[mu] = new ACECTildeBasisFunction[tot_size];
- }
-
-
- for (SPECIES_TYPE mu = 0; mu < this->nelements; ++mu) {
- size_t ind_rank1 = 0;
- size_t ind = 0;
-
- for (auto &func: mu0_ctilde_basis_vector[mu]) {
- if (func.rank == 1) { //r=0, rank=1
- basis_rank1[mu][ind_rank1] = func;
- ind_rank1 += 1;
- } else { //r>0, rank>1
- basis[mu][ind] = func;
- ind += 1;
- }
- }
-
- }
-}
-
-
-void ACECTildeBasisSet::_load_radial_ACERadial(FILE *fptr,
- const string filename,
- const string radbasename) {
- int res;
- char buffer[1024], buffer2[1024];
-
- res = fscanf(fptr, " nradbase=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "nradbase", "nradbase=[number]");
-// throw invalid_argument(("File '" + filename + "': couldn't read nradbase").c_str());
- nradbase = stoi_err(buffer, filename, "nradbase", "nradbase=[number]");
-
- res = fscanf(fptr, " nradmax=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "nradmax", "nradmax=[number]");
- nradmax = stoi_err(buffer, filename, "nradmax", "nradmax=[number]");
-
- res = fscanf(fptr, " cutoffmax=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "cutoffmax", "cutoffmax=[number]");
- cutoffmax = stod_err(buffer, filename, "cutoffmax", "cutoffmax=[number]");
-
-
- res = fscanf(fptr, " deltaSplineBins=");
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "deltaSplineBins", "deltaSplineBins=[spline density, Angstroms]");
-// throw invalid_argument(("File '" + filename + "': couldn't read ntot").c_str());
- deltaSplineBins = stod_err(buffer, filename, "deltaSplineBins", "deltaSplineBins=[spline density, Angstroms]");
-
-
- if (radial_functions == nullptr)
- radial_functions = new ACERadialFunctions(nradbase, lmax, nradmax,
- deltaSplineBins,
- nelements,
- cutoffmax, radbasename);
- else
- radial_functions->init(nradbase, lmax, nradmax,
- deltaSplineBins,
- nelements,
- cutoffmax, radbasename);
-
-
- //hard-core repulsion
- res = fscanf(fptr, " core repulsion parameters:");
- if (res != 0)
- throw_error(filename, "core repulsion parameters", "core repulsion parameters: [prehc lambdahc] ...");
-// throw invalid_argument(("File '" + filename + "': couldn't read core repulsion parameters").c_str());
-
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
- res = fscanf(fptr, "%s %s", buffer, buffer2);
- if (res != 2)
- throw_error(filename, "core repulsion parameters", "core repulsion parameters: [prehc lambdahc] ...");
- radial_functions->prehc(mu_i, mu_j) = stod_err(buffer, filename, "core repulsion parameters",
- "core repulsion parameters: [prehc lambdahc] ...");
- radial_functions->lambdahc(mu_i, mu_j) = stod_err(buffer2, filename, "core repulsion parameters",
- "core repulsion parameters: [prehc lambdahc] ...");
- }
-
-
-
- //read radial functions parameter
- res = fscanf(fptr, " radparameter=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "radparameter", "radparameter=[param_ele1] [param_ele2]");
- radial_functions->lambda(mu_i, mu_j) = stod_err(buffer, filename, "radparameter",
- "radparameter=[param_ele1] [param_ele2]");
- }
-
-
- res = fscanf(fptr, " cutoff=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "cutoff", "cutoff=[param_ele1] [param_ele2]");
- radial_functions->cut(mu_i, mu_j) = stod_err(buffer, filename, "cutoff",
- "cutoff=[param_ele1] [param_ele2]");
- }
-
-
- res = fscanf(fptr, " dcut=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
- res = fscanf(fptr, " %s", buffer);
- if (res != 1)
- throw_error(filename, "dcut", "dcut=[param_ele1] [param_ele2]");
- radial_functions->dcut(mu_i, mu_j) = stod_err(buffer, filename, "dcut", "dcut=[param_ele1] [param_ele2]");
- }
-
-
- res = fscanf(fptr, " crad=");
- for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i)
- for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j)
- for (NS_TYPE k = 0; k < nradbase; k++)
- for (NS_TYPE n = 0; n < nradmax; n++)
- for (LS_TYPE l = 0; l <= lmax; l++) {
- res = fscanf(fptr, "%s", buffer);
- if (res != 1)
- throw_error(filename, "crad", "crad=[crad_]...[crad_knl]: nradbase*nrad*(l+1) times");
- radial_functions->crad(mu_i, mu_j, n, l, k) = stod_err(buffer, filename, "crad",
- "crad=[crad_]...[crad_knl]: nradbase*nrad*(l+1) times");
- }
-}
-
-void ACECTildeBasisSet::_load_radial_SHIPsBasic(FILE *fptr,
- const string filename,
- const string radbasename) {
- // create a radial basis object, and read it from the file pointer
- SHIPsRadialFunctions *ships_radial_functions = new SHIPsRadialFunctions();
- ships_radial_functions->fread(fptr);
-
- //mimic ships_radial_functions to ACERadialFunctions
- ships_radial_functions->nradial = ships_radial_functions->get_maxn();
- ships_radial_functions->nradbase = ships_radial_functions->get_maxn();
-
- nradbase = ships_radial_functions->get_maxn();
- nradmax = ships_radial_functions->get_maxn();
- cutoffmax = ships_radial_functions->get_rcut();
- deltaSplineBins = 0.001;
-
- ships_radial_functions->init(nradbase, lmax, nradmax,
- deltaSplineBins,
- nelements,
- cutoffmax, radbasename);
-
-
- if (radial_functions) delete radial_functions;
- radial_functions = ships_radial_functions;
- radial_functions->prehc.fill(0);
- radial_functions->lambdahc.fill(1);
- radial_functions->lambda.fill(0);
-
-
- radial_functions->cut.fill(ships_radial_functions->get_rcut());
- radial_functions->dcut.fill(0);
-
- radial_functions->crad.fill(0);
-
-}
\ No newline at end of file
diff --git a/lib/pace/ace_c_basis.h b/lib/pace/ace_c_basis.h
deleted file mode 100644
index ec6b9e8afc..0000000000
--- a/lib/pace/ace_c_basis.h
+++ /dev/null
@@ -1,155 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Yury Lysogorskiy on 1.04.20.
-
-#ifndef ACE_C_BASIS_H
-#define ACE_C_BASIS_H
-
-#include "ace_flatten_basis.h"
-
-typedef vector> C_tilde_full_basis_vector2d;
-
-/**
- * ACE basis set of C-tilde basis functions
- */
-class ACECTildeBasisSet : public ACEFlattenBasisSet {
-public:
-
- ACECTildeBasisFunction **basis_rank1 = nullptr; ///< two-dimensional array of first-rank basis function with indices: [species index][func index]
- ACECTildeBasisFunction **basis = nullptr; ///< two-dimensional array of higher rank basis function with indices: [species index][func index]
-
- DOUBLE_TYPE *full_c_tildes_rank1 = nullptr; ///< C_tilde coefficients contiguous package, size: coeff_array_total_size_rank1
- DOUBLE_TYPE *full_c_tildes = nullptr; ///< C_tilde coefficients contiguous package, size: coeff_array_total_size
-
- //TODO: remove?
- SHORT_INT_TYPE num_ctilde_max = 0;
-
-
- /**
- * Default constructor
- */
- ACECTildeBasisSet() = default;
-
- /**
- * Constructor from .ace file
- */
- ACECTildeBasisSet(const string filename);
-
- /**
- * Copy constructor (see. Rule of Three)
- * @param other
- */
- ACECTildeBasisSet(const ACECTildeBasisSet &other);
-
- /**
- * operator= (see. Rule of Three)
- * @param other
- * @return
- */
- ACECTildeBasisSet &operator=(const ACECTildeBasisSet &other);
-
- /**
- * Destructor (see. Rule of Three)
- */
- ~ACECTildeBasisSet() override;
-
- /**
- * Cleaning dynamic memory of the class (see. Rule of Three)
- */
- void _clean() override;
-
- /**
- * Copying and cleaning dynamic memory of the class (see. Rule of Three)
- * @param src
- */
- void _copy_dynamic_memory(const ACECTildeBasisSet &src);
-
- /**
- * Copying scalar variables
- * @param src
- */
- void _copy_scalar_memory(const ACECTildeBasisSet &src);
-
- /**
- * Clean contiguous arrays (full_c_tildes_rank1, full_c_tildes) and those of base class
- */
- void _clean_contiguous_arrays() override ;
-
- /**
- * Save potential to .ace file
- * @param filename .ace file name
- */
- void save(const string &filename) override;
-
- /**
- * Load potential from .ace
- * @param filename .ace file name
- */
- void load(const string filename) override;
-
- /**
- * Load the ACE type radial basis
- */
- void _load_radial_ACERadial(FILE *fptr,
- const string filename,
- const string radbasename);
-
- void _load_radial_SHIPsBasic(FILE * fptr,
- const string filename,
- const string radbasename );
-
- /**
- * Re-pack the constituent dynamic arrays of all basis functions in contiguous arrays
- */
- void pack_flatten_basis() override;
-
- /**
- * Computes flatten array sizes
- * @param basis_rank1 two-dimensional array of first-rank ACECTildeBasisFunctions
- * @param basis two-dimensional array of higher-rank ACECTildeBasisFunctions
- */
- void compute_array_sizes(ACECTildeBasisFunction** basis_rank1, ACECTildeBasisFunction** basis);
-
- /**
- * Clean basis arrays 'basis_rank1' and 'basis'
- */
- void _clean_basis_arrays();
-
- /**
- * Pack two-dimensional vector of ACECTildeBasisFunction into 1D dynami array with all basis functions
- * @param mu0_ctilde_basis_vector vector>
- */
- void flatten_basis(C_tilde_full_basis_vector2d& mu0_ctilde_basis_vector);
-
- /**
- * Empty stub implementation
- */
- void flatten_basis() override{};
-};
-
-#endif //ACE_C_BASIS_H
diff --git a/lib/pace/ace_c_basisfunction.h b/lib/pace/ace_c_basisfunction.h
deleted file mode 100644
index 4e2795590e..0000000000
--- a/lib/pace/ace_c_basisfunction.h
+++ /dev/null
@@ -1,251 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Yury Lysogorskiy on 26.02.20.
-
-#ifndef ACE_C_BASISFUNCTION_H
-#define ACE_C_BASISFUNCTION_H
-
-#include
-#include
-#include
-#include
-
-#include "ace_types.h"
-
-//macros for copying the member-array from "other" object for C-tilde and B-basis
-#define basis_mem_copy(other, array, size, type) if(other.array) { \
- if(!is_proxy) delete[] array;\
- array = new type[(size)];\
- is_proxy = false;\
- memcpy(array, other.array, (size) * sizeof(type));\
-}
-
-/**
- * Abstract basis function, that stores general quantities
- */
-struct ACEAbstractBasisFunction {
- /**
- * flattened array of computed combinations of (m1, m2, ..., m_rank)
- * which have non-zero general Clebsch-Gordan coefficient:
- * \f$ \mathbf{m}_1, \dots, \mathbf{m}_\mathrm{num ms combs}\f$ =
- * \f$ (m_1, m_2, \dots, m_{rank})_1, \dots, (m_1, m_2, \dots, m_{rank})_{\mathrm{num ms combs}} \f$,
- * size = num_ms_combs * rank,
- * effective shape: [num_ms_combs][rank]
- */
- MS_TYPE *ms_combs = nullptr;
-
- /**
- * species types of neighbours atoms \f$ \mathbf{\mu} = (\mu_1, \mu_2, \dots, \mu_{rank}) \f$,
- * should be lexicographically sorted,
- * size = rank,
- * effective shape: [rank]
- */
- SPECIES_TYPE *mus = nullptr;
-
- /**
- * indices for radial part \f$ \mathbf{n} = (n_1, n_2, \dots, n_{rank}) \f$,
- * should be lexicographically sorted,
- * size = rank,
- * effective shape: [rank]
- */
- NS_TYPE *ns = nullptr;
-
-
- /**
- * indices for radial part \f$ \mathbf{l} = (l_1, l_2, \dots, l_{rank}) \f$,
- * should be lexicographically sorted,
- * size = rank,
- * effective shape: [rank]
- */
- LS_TYPE *ls = nullptr;
-
- SHORT_INT_TYPE num_ms_combs = 0; ///< number of different ms-combinations
-
- RANK_TYPE rank = 0; ///< number of atomic base functions "A"s in basis function product B
-
- DENSITY_TYPE ndensity = 0; ///< number of densities
-
- SPECIES_TYPE mu0 = 0; ///< species type of central atom
-
- /**
- * whether ms array contains only "non-negative" ms-combinations.
- * positive ms-combination is when the first non-zero m is positive (0 1 -1)
- * negative ms-combination is when the first non-zero m is negative (0 -1 1)
- */
- bool is_half_ms_basis = false;
-
- /*
- * flag, whether object is "owner" (i.e. responsible for memory cleaning) of
- * the ms, ns, ls, mus and other dynamically allocated arrases or just a proxy for them
- */
- bool is_proxy = false;
-
- /**
- * Copying static and dynamic memory variables from another ACEAbstractBasisFunction.
- * Always copy the dynamic memory, even if the source is a proxy object
- * @param other
- */
- virtual void _copy_from(const ACEAbstractBasisFunction &other) {
- rank = other.rank;
- ndensity = other.ndensity;
- mu0 = other.mu0;
- num_ms_combs = other.num_ms_combs;
- is_half_ms_basis = other.is_half_ms_basis;
- is_proxy = false;
-
- basis_mem_copy(other, mus, rank, SPECIES_TYPE)
- basis_mem_copy(other, ns, rank, NS_TYPE)
- basis_mem_copy(other, ls, rank, LS_TYPE)
- basis_mem_copy(other, ms_combs, num_ms_combs * rank, MS_TYPE)
- }
-
- /**
- * Clean the dynamically allocated memory if object is responsible for it
- */
- virtual void _clean() {
- //release memory if the structure is not proxy
- if (!is_proxy) {
- delete[] mus;
- delete[] ns;
- delete[] ls;
- delete[] ms_combs;
- }
-
- mus = nullptr;
- ns = nullptr;
- ls = nullptr;
- ms_combs = nullptr;
- }
-
-};
-
-/**
- * Representation of C-tilde basis function, i.e. the function that is summed up over a group of B-functions
- * that differs only by intermediate coupling orbital moment \f$ L \f$ and coefficients.
- */
-struct ACECTildeBasisFunction : public ACEAbstractBasisFunction {
-
- /**
- * c_tilde coefficients for all densities,
- * size = num_ms_combs*ndensity,
- * effective shape [num_ms_combs][ndensity]
- */
- DOUBLE_TYPE *ctildes = nullptr;
-
- /*
- * Default constructor
- */
- ACECTildeBasisFunction() = default;
-
- // Because the ACECTildeBasisFunction contains dynamically allocated arrays, the Rule of Three should be
- // fulfilled, i.e. copy constructor (copy the dynamic arrays), operator= (release previous arrays and
- // copy the new dynamic arrays) and destructor (release all dynamically allocated memory)
-
- /**
- * Copy constructor, to fulfill the Rule of Three.
- * Always copy the dynamic memory, even if the source is a proxy object.
- */
- ACECTildeBasisFunction(const ACECTildeBasisFunction &other) {
- _copy_from(other);
- }
-
- /*
- * operator=, to fulfill the Rule of Three.
- * Always copy the dynamic memory, even if the source is a proxy object
- */
- ACECTildeBasisFunction &operator=(const ACECTildeBasisFunction &other) {
- _clean();
- _copy_from(other);
- return *this;
- }
-
- /*
- * Destructor
- */
- ~ACECTildeBasisFunction() {
- _clean();
- }
-
- /**
- * Copy from another object, always copy the memory, even if the source is a proxy object
- * @param other
- */
- void _copy_from(const ACECTildeBasisFunction &other) {
- ACEAbstractBasisFunction::_copy_from(other);
- is_proxy = false;
- basis_mem_copy(other, ctildes, num_ms_combs * ndensity, DOUBLE_TYPE)
- }
-
- /**
- * Clean the dynamically allocated memory if object is responsible for it
- */
- void _clean() override {
- ACEAbstractBasisFunction::_clean();
- //release memory if the structure is not proxy
- if (!is_proxy) {
- delete[] ctildes;
- }
- ctildes = nullptr;
- }
-
- /**
- * Print the information about basis function to cout, in order to ease the output redirection.
- */
- void print() const {
- using namespace std;
- cout << "ACECTildeBasisFunction: ndensity= " << (int) this->ndensity << ", mu0 = " << (int) this->mu0 << " ";
- cout << " mus=(";
- for (RANK_TYPE r = 0; r < this->rank; r++)
- cout << (int) this->mus[r] << " ";
- cout << "), ns=(";
- for (RANK_TYPE r = 0; r < this->rank; r++)
- cout << this->ns[r] << " ";
- cout << "), ls=(";
- for (RANK_TYPE r = 0; r < this->rank; r++)
- cout << this->ls[r] << " ";
-
- cout << "), " << this->num_ms_combs << " m_s combinations: {" << endl;
-
- for (int i = 0; i < this->num_ms_combs; i++) {
- cout << "\t< ";
- for (RANK_TYPE r = 0; r < this->rank; r++)
- cout << this->ms_combs[i * this->rank + r] << " ";
- cout << " >: c_tilde=";
- for (DENSITY_TYPE p = 0; p < this->ndensity; ++p)
- cout << " " << this->ctildes[i * this->ndensity + p] << " ";
- cout << endl;
- }
- if (this->is_proxy)
- cout << "proxy ";
- if (this->is_half_ms_basis)
- cout << "half_ms_basis";
- cout << "}" << endl;
- }
-};
-
-#endif //ACE_C_BASISFUNCTION_H
diff --git a/lib/pace/ace_complex.h b/lib/pace/ace_complex.h
deleted file mode 100644
index 5fdb4b5b93..0000000000
--- a/lib/pace/ace_complex.h
+++ /dev/null
@@ -1,266 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Yury Lysogorskiy on 26.02.20.
-
-#ifndef ACE_COMPLEX_H
-#define ACE_COMPLEX_H
-
-
-/**
-A custom data structure for complex numbers and overloaded operations with them.
-*/
-struct ACEComplex {
-public:
- /**
- Double, real part of the complex number
- */
- DOUBLE_TYPE real;
- /**
- Double, imaginary part of the complex number
- */
- DOUBLE_TYPE img;
-
- ACEComplex &operator=(const ACEComplex &rhs) = default;
-
- ACEComplex &operator=(const DOUBLE_TYPE &rhs) {
- this->real = rhs;
- this->img = 0.;
- return *this;
- }
-
- /**
- Overloading of arithmetical operator += ACEComplex
- */
- ACEComplex &operator+=(const ACEComplex &rhs) {
- this->real += rhs.real;
- this->img += rhs.img;
- return *this; // return the result by reference
- }
-
- /**
- Overloading of arithmetical operator += DOUBLE_TYPE
- */
- ACEComplex &operator+=(const DOUBLE_TYPE &rhs) {
- this->real += rhs;
- return *this; // return the result by reference
- }
-
- /**
- Overloading of arithmetical operator *= DOUBLE_TYPE
- */
- ACEComplex &operator*=(const DOUBLE_TYPE &rhs) {
- this->real *= rhs;
- this->img *= rhs;
- return *this; // return the result by reference
- }
-
- /**
- Overloading of arithmetical operator *= ACEComplex
- */
- ACEComplex &operator*=(const ACEComplex &rhs) {
- DOUBLE_TYPE old_real = this->real;
- this->real = old_real * rhs.real - this->img * rhs.img;
- this->img = old_real * rhs.img + this->img * rhs.real;
- return *this; // return the result by reference
- }
-
- /**
- Overloading of arithmetical operator * ACEComplex
- */
- ACEComplex operator*(const ACEComplex &cm2) const {
- ACEComplex res{real * cm2.real - img * cm2.img,
- real * cm2.img + img * cm2.real};
- return res;
- }
-
- /*
- * Return complex conjugated copy of itself
- */
- ACEComplex conjugated() const {
- ACEComplex res{real, -img};
- return res;
- }
-
- /*
- * Complex conjugate itself inplace
- */
- void conjugate() {
- img = -img;
- }
-
- /*
- * Multiplication by ACEComplex and return real-part only
- */
- DOUBLE_TYPE real_part_product(const ACEComplex &cm2) const {
- return real * cm2.real - img * cm2.img;
- }
-
- /*
- * Multiplication by DOUBLE_TYPE and return real-part only
- */
- DOUBLE_TYPE real_part_product(const DOUBLE_TYPE &cm2) const {
- return real * cm2;
- }
-
- /*
- * Overloading of arithmetical operator * by DOUBLE_TYPE
- */
- ACEComplex operator*(const DOUBLE_TYPE &cm2) const {
- ACEComplex res{real * cm2,
- img * cm2};
- return res;
- }
-
- /*
- * Overloading of arithmetical operator + ACEComplex
- */
- ACEComplex operator+(const ACEComplex &cm2) const {
- ACEComplex res{real + cm2.real,
- img + cm2.img};
- return res;
- }
-
- /*
- * Overloading of arithmetical operator + with DOUBLE_TYPE
- */
- ACEComplex operator+(const DOUBLE_TYPE &cm2) const {
- ACEComplex res{real + cm2, img};
- return res;
- }
-
- /*
- * Overloading of arithmetical operator == ACEComplex
- */
- bool operator==(const ACEComplex &c2) const {
- return (real == c2.real) && (img == c2.img);
- }
-
- /*
- * Overloading of arithmetical operator == DOUBLE_TYPE
- */
- bool operator==(const DOUBLE_TYPE &d2) const {
- return (real == d2) && (img == 0.);
- }
-
- /*
- * Overloading of arithmetical operator != ACEComplex
- */
- bool operator!=(const ACEComplex &c2) const {
- return (real != c2.real) || (img != c2.img);
- }
-
- /*
- * Overloading of arithmetical operator != DOUBLE_TYPE
- */
- bool operator!=(const DOUBLE_TYPE &d2) const {
- return (real != d2) || (img != 0.);
- }
-
-};
-
-/*
- * double * complex for commutativity with complex * double
- */
-inline ACEComplex operator*(const DOUBLE_TYPE &real, const ACEComplex &cm) {
- return cm * real;
-}
-
-/*
- * double + complex for commutativity with complex + double
- */
-inline ACEComplex operator+(const DOUBLE_TYPE &real, const ACEComplex &cm) {
- return cm + real;
-}
-
-/**
-A structure to store the derivative of \f$ Y_{lm} \f$
-*/
-struct ACEDYcomponent {
-public:
- /**
- complex, contains the three components of derivative of Ylm,
- \f$ \frac{dY_{lm}}{dx}, \frac{dY_{lm}}{dy} and \frac{dY_{lm}}{dz}\f$
- */
- ACEComplex a[3];
-
- /*
- * Overloading of arithmetical operator*= DOUBLE_TYPE
- */
- ACEDYcomponent &operator*=(const DOUBLE_TYPE &rhs) {
- this->a[0] *= rhs;
- this->a[1] *= rhs;
- this->a[2] *= rhs;
- return *this;
- }
-
- /*
- * Overloading of arithmetical operator * ACEComplex
- */
- ACEDYcomponent operator*(const ACEComplex &rhs) const {
- ACEDYcomponent res;
- res.a[0] = this->a[0] * rhs;
- res.a[1] = this->a[1] * rhs;
- res.a[2] = this->a[2] * rhs;
- return res;
- }
-
- /*
- * Overloading of arithmetical operator * DOUBLE_TYPE
- */
- ACEDYcomponent operator*(const DOUBLE_TYPE &rhs) const {
- ACEDYcomponent res;
- res.a[0] = this->a[0] * rhs;
- res.a[1] = this->a[1] * rhs;
- res.a[2] = this->a[2] * rhs;
- return res;
- }
-
- /*
- * Return conjugated copy of itself
- */
- ACEDYcomponent conjugated() const {
- ACEDYcomponent res;
- res.a[0] = this->a[0].conjugated();
- res.a[1] = this->a[1].conjugated();
- res.a[2] = this->a[2].conjugated();
- return res;
- }
-
- /*
- * Conjugated itself in-place
- */
- void conjugate() {
- this->a[0].conjugate();
- this->a[1].conjugate();
- this->a[2].conjugate();
- }
-
-};
-
-
-#endif //ACE_COMPLEX_H
diff --git a/lib/pace/ace_contigous_array.h b/lib/pace/ace_contigous_array.h
deleted file mode 100644
index f008fa203f..0000000000
--- a/lib/pace/ace_contigous_array.h
+++ /dev/null
@@ -1,249 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-
-// Created by Yury Lysogorskiy on 11.01.20.
-#ifndef ACE_CONTIGUOUSARRAYND_H
-#define ACE_CONTIGUOUSARRAYND_H
-
-#include
-
-#include "ace_types.h"
-
-using namespace std;
-
-/**
- * Common predecessor class to represent multidimensional array of type T
- * and store it in memory contiguous form
- *
- * @tparam T data type
- */
-template
-class ContiguousArrayND {
-protected:
- T *data = nullptr; ///< pointer to contiguous data
- size_t size = 0; ///< total array size
- string array_name = "Array"; /// 0) {
- data = new T[size];
- for (size_t ind = 0; ind < size; ind++)
- data[ind] = other.data[ind];
- }
- } else { //is proxy, then copy the pointer
- data = other.data;
- }
- }
-
- /**
- * Overload operator=
- * @param other another ContiguousArrayND
- * @return itself
- */
-
- ContiguousArrayND &operator=(const ContiguousArrayND &other) {
-#ifdef MULTIARRAY_LIFE_CYCLE
- cout< 0) {
-
- if(data!=nullptr) delete[] data;
- data = new T[size];
-
- for (size_t ind = 0; ind < size; ind++)
- data[ind] = other.data[ind];
- }
- } else { //is proxy, then copy the pointer
- data = other.data;
- }
- }
- return *this;
- }
-
-
- //TODO: make destructor virtual, check the destructors in inherited classes
-
- /**
- * Destructor
- */
- ~ContiguousArrayND() {
-#ifdef MULTIARRAY_LIFE_CYCLE
- cout<array_name = name;
- }
-
- /**
- * Get total number of elements in array (its size)
- * @return array size
- */
- size_t get_size() const {
- return size;
- }
-
- /**
- * Fill array with value
- * @param value value to fill
- */
- void fill(T value) {
- for (size_t ind = 0; ind < size; ind++)
- data[ind] = value;
- }
-
- /**
- * Get array data at absolute index ind for reading
- * @param ind absolute index
- * @return array value
- */
- inline const T &get_data(size_t ind) const {
-#ifdef MULTIARRAY_INDICES_CHECK
- if ((ind < 0) | (ind >= size)) {
- printf("%s: get_data ind=%d out of range (0, %d)\n", array_name, ind, size);
- exit(EXIT_FAILURE);
- }
-#endif
- return data[ind];
- }
-
- /**
- * Get array data at absolute index ind for writing
- * @param ind absolute index
- * @return array value
- */
- inline T &get_data(size_t ind) {
-#ifdef MULTIARRAY_INDICES_CHECK
- if ((ind < 0) | (ind >= size)) {
- printf("%s: get_data ind=%ld out of range (0, %ld)\n", array_name.c_str(), ind, size);
- exit(EXIT_FAILURE);
- }
-#endif
- return data[ind];
- }
-
- /**
- * Get array data pointer
- * @return data array pointer
- */
- inline T* get_data() const {
- return data;
- }
-
- /**
- * Overload comparison operator==
- * Compare the total size and array values elementwise.
- *
- * @param other another array
- * @return
- */
- bool operator==(const ContiguousArrayND &other) const {
- if (this->size != other.size)
- return false;
-
- for (size_t i = 0; i < this->size; ++i)
- if (this->data[i] != other.data[i])
- return false;
-
- return true;
- }
-
-
- /**
- * Convert to flatten vector container
- * @return vector container
- */
- vector to_flatten_vector() const {
- vector res;
-
- res.resize(size);
- size_t vec_ind = 0;
-
- for (int vec_ind = 0; vec_ind < size; vec_ind++)
- res.at(vec_ind) = data[vec_ind];
-
- return res;
- } // end to_flatten_vector()
-
-
- /**
- * Set values from flatten vector container
- * @param vec container
- */
- void set_flatten_vector(const vector &vec) {
- if (vec.size() != size)
- throw std::invalid_argument("Flatten vector size is not consistent with expected size");
- for (size_t i = 0; i < size; i++) {
- data[i] = vec[i];
- }
- }
-
- bool is_proxy(){
- return is_proxy_;
- }
-
-};
-
-
-#endif //ACE_CONTIGUOUSARRAYND_H
\ No newline at end of file
diff --git a/lib/pace/ace_evaluator.cpp b/lib/pace/ace_evaluator.cpp
deleted file mode 100644
index 987f4502bf..0000000000
--- a/lib/pace/ace_evaluator.cpp
+++ /dev/null
@@ -1,660 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-// Created by Yury Lysogorskiy on 31.01.20.
-
-#include "ace_evaluator.h"
-
-#include "ace_abstract_basis.h"
-#include "ace_types.h"
-
-void ACEEvaluator::init(ACEAbstractBasisSet *basis_set) {
- A.init(basis_set->nelements, basis_set->nradmax + 1, basis_set->lmax + 1, "A");
- A_rank1.init(basis_set->nelements, basis_set->nradbase, "A_rank1");
-
- rhos.init(basis_set->ndensitymax + 1, "rhos"); // +1 density for core repulsion
- dF_drho.init(basis_set->ndensitymax + 1, "dF_drho"); // +1 density for core repulsion
-}
-
-void ACEEvaluator::init_timers() {
- loop_over_neighbour_timer.init();
- forces_calc_loop_timer.init();
- forces_calc_neighbour_timer.init();
- energy_calc_timer.init();
- per_atom_calc_timer.init();
- total_time_calc_timer.init();
-}
-
-//================================================================================================================
-
-void ACECTildeEvaluator::set_basis(ACECTildeBasisSet &bas) {
- basis_set = &bas;
- init(basis_set);
-}
-
-void ACECTildeEvaluator::init(ACECTildeBasisSet *basis_set) {
-
- ACEEvaluator::init(basis_set);
-
-
- weights.init(basis_set->nelements, basis_set->nradmax + 1, basis_set->lmax + 1,
- "weights");
-
- weights_rank1.init(basis_set->nelements, basis_set->nradbase, "weights_rank1");
-
-
- DG_cache.init(1, basis_set->nradbase, "DG_cache");
- DG_cache.fill(0);
-
- R_cache.init(1, basis_set->nradmax, basis_set->lmax + 1, "R_cache");
- R_cache.fill(0);
-
- DR_cache.init(1, basis_set->nradmax, basis_set->lmax + 1, "DR_cache");
- DR_cache.fill(0);
-
- Y_cache.init(1, basis_set->lmax + 1, "Y_cache");
- Y_cache.fill({0, 0});
-
- DY_cache.init(1, basis_set->lmax + 1, "dY_dense_cache");
- DY_cache.fill({0.});
-
- //hard-core repulsion
- DCR_cache.init(1, "DCR_cache");
- DCR_cache.fill(0);
- dB_flatten.init(basis_set->max_dB_array_size, "dB_flatten");
-
-
-}
-
-void ACECTildeEvaluator::resize_neighbours_cache(int max_jnum) {
- if(basis_set== nullptr) {
- throw std::invalid_argument("ACECTildeEvaluator: basis set is not assigned");
- }
- if (R_cache.get_dim(0) < max_jnum) {
-
- //TODO: implement grow
- R_cache.resize(max_jnum, basis_set->nradmax, basis_set->lmax + 1);
- R_cache.fill(0);
-
- DR_cache.resize(max_jnum, basis_set->nradmax, basis_set->lmax + 1);
- DR_cache.fill(0);
-
- DG_cache.resize(max_jnum, basis_set->nradbase);
- DG_cache.fill(0);
-
- Y_cache.resize(max_jnum, basis_set->lmax + 1);
- Y_cache.fill({0, 0});
-
- DY_cache.resize(max_jnum, basis_set->lmax + 1);
- DY_cache.fill({0});
-
- //hard-core repulsion
- DCR_cache.init(max_jnum, "DCR_cache");
- DCR_cache.fill(0);
- }
-}
-
-
-
-// double** r - atomic coordinates of atom I
-// int* types - atomic types if atom I
-// int **firstneigh - ptr to 1st J int value of each I atom. Usage: jlist = firstneigh[i];
-// Usage: j = jlist_of_i[jj];
-// jnum - number of J neighbors for each I atom. jnum = numneigh[i];
-
-void
-ACECTildeEvaluator::compute_atom(int i, DOUBLE_TYPE **x, const SPECIES_TYPE *type, const int jnum, const int *jlist) {
- if(basis_set== nullptr) {
- throw std::invalid_argument("ACECTildeEvaluator: basis set is not assigned");
- }
- per_atom_calc_timer.start();
-#ifdef PRINT_MAIN_STEPS
- printf("\n ATOM: ind = %d r_norm=(%f, %f, %f)\n",i, x[i][0], x[i][1], x[i][2]);
-#endif
- DOUBLE_TYPE evdwl = 0, evdwl_cut = 0, rho_core = 0;
- DOUBLE_TYPE r_norm;
- DOUBLE_TYPE xn, yn, zn, r_xyz;
- DOUBLE_TYPE R, GR, DGR, R_over_r, DR, DCR;
- DOUBLE_TYPE *r_hat;
-
- SPECIES_TYPE mu_j;
- RANK_TYPE r, rank, t;
- NS_TYPE n;
- LS_TYPE l;
- MS_TYPE m, m_t;
-
- SPECIES_TYPE *mus;
- NS_TYPE *ns;
- LS_TYPE *ls;
- MS_TYPE *ms;
-
- int j, jj, func_ind, ms_ind;
- SHORT_INT_TYPE factor;
-
- ACEComplex Y{0}, Y_DR{0.};
- ACEComplex B{0.};
- ACEComplex dB{0};
- ACEComplex A_cache[basis_set->rankmax];
-
- dB_flatten.fill({0.});
-
- ACEDYcomponent grad_phi_nlm{0}, DY{0.};
-
- //size is +1 of max to avoid out-of-boundary array access in double-triangular scheme
- ACEComplex A_forward_prod[basis_set->rankmax + 1];
- ACEComplex A_backward_prod[basis_set->rankmax + 1];
-
- DOUBLE_TYPE inv_r_norm;
- DOUBLE_TYPE r_norms[jnum];
- DOUBLE_TYPE inv_r_norms[jnum];
- DOUBLE_TYPE rhats[jnum][3];//normalized vector
- SPECIES_TYPE elements[jnum];
- const DOUBLE_TYPE xtmp = x[i][0];
- const DOUBLE_TYPE ytmp = x[i][1];
- const DOUBLE_TYPE ztmp = x[i][2];
- DOUBLE_TYPE f_ji[3];
-
- bool is_element_mapping = element_type_mapping.get_size() > 0;
- SPECIES_TYPE mu_i;
- if (is_element_mapping)
- mu_i = element_type_mapping(type[i]);
- else
- mu_i = type[i];
-
- const SHORT_INT_TYPE total_basis_size_rank1 = basis_set->total_basis_size_rank1[mu_i];
- const SHORT_INT_TYPE total_basis_size = basis_set->total_basis_size[mu_i];
-
- ACECTildeBasisFunction *basis_rank1 = basis_set->basis_rank1[mu_i];
- ACECTildeBasisFunction *basis = basis_set->basis[mu_i];
-
- DOUBLE_TYPE rho_cut, drho_cut, fcut, dfcut;
- DOUBLE_TYPE dF_drho_core;
-
- //TODO: lmax -> lmaxi (get per-species type)
- const LS_TYPE lmaxi = basis_set->lmax;
-
- //TODO: nradmax -> nradiali (get per-species type)
- const NS_TYPE nradiali = basis_set->nradmax;
-
- //TODO: nradbase -> nradbasei (get per-species type)
- const NS_TYPE nradbasei = basis_set->nradbase;
-
- //TODO: get per-species type number of densities
- const DENSITY_TYPE ndensity= basis_set->ndensitymax;
-
- neighbours_forces.resize(jnum, 3);
- neighbours_forces.fill(0);
-
- //TODO: shift nullifications to place where arrays are used
- weights.fill({0});
- weights_rank1.fill(0);
- A.fill({0});
- A_rank1.fill(0);
- rhos.fill(0);
- dF_drho.fill(0);
-
-#ifdef EXTRA_C_PROJECTIONS
- basis_projections_rank1.init(total_basis_size_rank1, ndensity, "c_projections_rank1");
- basis_projections.init(total_basis_size, ndensity, "c_projections");
-#endif
-
- //proxy references to spherical harmonics and radial functions arrays
- const Array2DLM &ylm = basis_set->spherical_harmonics.ylm;
- const Array2DLM &dylm = basis_set->spherical_harmonics.dylm;
-
- const Array2D &fr = basis_set->radial_functions->fr;
- const Array2D &dfr = basis_set->radial_functions->dfr;
-
- const Array1D &gr = basis_set->radial_functions->gr;
- const Array1D &dgr = basis_set->radial_functions->dgr;
-
- loop_over_neighbour_timer.start();
-
- int jj_actual = 0;
- SPECIES_TYPE type_j = 0;
- int neighbour_index_mapping[jnum]; // jj_actual -> jj
- //loop over neighbours, compute distance, consider only atoms within with rradial_functions->cut(mu_i, mu_j);
- r_xyz = sqrt(xn * xn + yn * yn + zn * zn);
-
- if (r_xyz >= current_cutoff)
- continue;
-
- inv_r_norm = 1 / r_xyz;
-
- r_norms[jj_actual] = r_xyz;
- inv_r_norms[jj_actual] = inv_r_norm;
- rhats[jj_actual][0] = xn * inv_r_norm;
- rhats[jj_actual][1] = yn * inv_r_norm;
- rhats[jj_actual][2] = zn * inv_r_norm;
- elements[jj_actual] = mu_j;
- neighbour_index_mapping[jj_actual] = jj;
- jj_actual++;
- }
-
- int jnum_actual = jj_actual;
-
- //ALGORITHM 1: Atomic base A
- for (jj = 0; jj < jnum_actual; ++jj) {
- r_norm = r_norms[jj];
- mu_j = elements[jj];
- r_hat = rhats[jj];
-
- //proxies
- Array2DLM &Y_jj = Y_cache(jj);
- Array2DLM &DY_jj = DY_cache(jj);
-
-
- basis_set->radial_functions->evaluate(r_norm, basis_set->nradbase, nradiali, mu_i, mu_j);
- basis_set->spherical_harmonics.compute_ylm(r_hat[0], r_hat[1], r_hat[2], lmaxi);
- //loop for computing A's
- //rank = 1
- for (n = 0; n < basis_set->nradbase; n++) {
- GR = gr(n);
-#ifdef DEBUG_ENERGY_CALCULATIONS
- printf("-neigh atom %d\n", jj);
- printf("gr(n=%d)(r=%f) = %f\n", n, r_norm, gr(n));
- printf("dgr(n=%d)(r=%f) = %f\n", n, r_norm, dgr(n));
-#endif
- DG_cache(jj, n) = dgr(n);
- A_rank1(mu_j, n) += GR * Y00;
- }
- //loop for computing A's
- // for rank > 1
- for (n = 0; n < nradiali; n++) {
- auto &A_lm = A(mu_j, n);
- for (l = 0; l <= lmaxi; l++) {
- R = fr(n, l);
-#ifdef DEBUG_ENERGY_CALCULATIONS
- printf("R(nl=%d,%d)(r=%f)=%f\n", n + 1, l, r_norm, R);
-#endif
-
- DR_cache(jj, n, l) = dfr(n, l);
- R_cache(jj, n, l) = R;
-
- for (m = 0; m <= l; m++) {
- Y = ylm(l, m);
-#ifdef DEBUG_ENERGY_CALCULATIONS
- printf("Y(lm=%d,%d)=(%f, %f)\n", l, m, Y.real, Y.img);
-#endif
- A_lm(l, m) += R * Y; //accumulation sum over neighbours
- Y_jj(l, m) = Y;
- DY_jj(l, m) = dylm(l, m);
- }
- }
- }
-
- //hard-core repulsion
- rho_core += basis_set->radial_functions->cr;
- DCR_cache(jj) = basis_set->radial_functions->dcr;
-
- } //end loop over neighbours
-
- //complex conjugate A's (for NEGATIVE (-m) terms)
- // for rank > 1
- for (mu_j = 0; mu_j < basis_set->nelements; mu_j++) {
- for (n = 0; n < nradiali; n++) {
- auto &A_lm = A(mu_j, n);
- for (l = 0; l <= lmaxi; l++) {
- //fill in -m part in the outer loop using the same m <-> -m symmetry as for Ylm
- for (m = 1; m <= l; m++) {
- factor = m % 2 == 0 ? 1 : -1;
- A_lm(l, -m) = A_lm(l, m).conjugated() * factor;
- }
- }
- }
- } //now A's are constructed
- loop_over_neighbour_timer.stop();
-
- // ==================== ENERGY ====================
-
- energy_calc_timer.start();
-#ifdef EXTRA_C_PROJECTIONS
- basis_projections_rank1.fill(0);
- basis_projections.fill(0);
-#endif
-
- //ALGORITHM 2: Basis functions B with iterative product and density rho(p) calculation
- //rank=1
- for (int func_rank1_ind = 0; func_rank1_ind < total_basis_size_rank1; ++func_rank1_ind) {
- ACECTildeBasisFunction *func = &basis_rank1[func_rank1_ind];
-// ndensity = func->ndensity;
-#ifdef PRINT_LOOPS_INDICES
- printf("Num density = %d r = 0\n",(int) ndensity );
- print_C_tilde_B_basis_function(*func);
-#endif
- double A_cur = A_rank1(func->mus[0], func->ns[0] - 1);
-#ifdef DEBUG_ENERGY_CALCULATIONS
- printf("A_r=1(x=%d, n=%d)=(%f)\n", func->mus[0], func->ns[0], A_cur);
- printf(" coeff[0] = %f\n", func->ctildes[0]);
-#endif
- for (DENSITY_TYPE p = 0; p < ndensity; ++p) {
- //for rank=1 (r=0) only 1 ms-combination exists (ms_ind=0), so index of func.ctildes is 0..ndensity-1
- rhos(p) += func->ctildes[p] * A_cur;
-#ifdef EXTRA_C_PROJECTIONS
- //aggregate C-projections separately
- basis_projections_rank1(func_rank1_ind, p)+= func->ctildes[p] * A_cur;
-#endif
- }
- } // end loop for rank=1
-
- //rank>1
- int func_ms_ind = 0;
- int func_ms_t_ind = 0;// index for dB
-
- for (func_ind = 0; func_ind < total_basis_size; ++func_ind) {
- ACECTildeBasisFunction *func = &basis[func_ind];
- //TODO: check if func->ctildes are zero, then skip
-// ndensity = func->ndensity;
- rank = func->rank;
- r = rank - 1;
-#ifdef PRINT_LOOPS_INDICES
- printf("Num density = %d r = %d\n",(int) ndensity, (int)r );
- print_C_tilde_B_basis_function(*func);
-#endif
- mus = func->mus;
- ns = func->ns;
- ls = func->ls;
-
- //loop over {ms} combinations in sum
- for (ms_ind = 0; ms_ind < func->num_ms_combs; ++ms_ind, ++func_ms_ind) {
- ms = &func->ms_combs[ms_ind * rank]; // current ms-combination (of length = rank)
-
- //loop over m, collect B = product of A with given ms
- A_forward_prod[0] = 1;
- A_backward_prod[r] = 1;
-
- //fill forward A-product triangle
- for (t = 0; t < rank; t++) {
- //TODO: optimize ns[t]-1 -> ns[t] during functions construction
- A_cache[t] = A(mus[t], ns[t] - 1, ls[t], ms[t]);
-#ifdef DEBUG_ENERGY_CALCULATIONS
- printf("A(x=%d, n=%d, l=%d, m=%d)=(%f,%f)\n", mus[t], ns[t], ls[t], ms[t], A_cache[t].real,
- A_cache[t].img);
-#endif
- A_forward_prod[t + 1] = A_forward_prod[t] * A_cache[t];
- }
-
- B = A_forward_prod[t];
-
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("B = (%f, %f)\n", (B).real, (B).img);
-#endif
- //fill backward A-product triangle
- for (t = r; t >= 1; t--) {
- A_backward_prod[t - 1] =
- A_backward_prod[t] * A_cache[t];
- }
-
- for (t = 0; t < rank; ++t, ++func_ms_t_ind) {
- dB = A_forward_prod[t] * A_backward_prod[t]; //dB - product of all A's except t-th
- dB_flatten(func_ms_t_ind) = dB;
-#ifdef DEBUG_FORCES_CALCULATIONS
- m_t = ms[t];
- printf("dB(n,l,m)(%d,%d,%d) = (%f, %f)\n", ns[t], ls[t], m_t, (dB).real, (dB).img);
-#endif
- }
-
- for (DENSITY_TYPE p = 0; p < ndensity; ++p) {
- //real-part only multiplication
- rhos(p) += B.real_part_product(func->ctildes[ms_ind * ndensity + p]);
-
-#ifdef EXTRA_C_PROJECTIONS
- //aggregate C-projections separately
- basis_projections(func_ind, p)+=B.real_part_product(func->ctildes[ms_ind * ndensity + p]);
-#endif
-
-#ifdef PRINT_INTERMEDIATE_VALUES
- printf("rhos(%d) += %f\n", p, B.real_part_product(func->ctildes[ms_ind * ndensity + p]));
- printf("Rho[i = %d][p = %d] = %f\n", i , p , rhos(p));
-#endif
- }
- }//end of loop over {ms} combinations in sum
- }// end loop for rank>1
-
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("rhos = ");
- for(DENSITY_TYPE p =0; prho_core_cutoffs(mu_i);
- drho_cut = basis_set->drho_core_cutoffs(mu_i);
-
- basis_set->inner_cutoff(rho_core, rho_cut, drho_cut, fcut, dfcut);
- basis_set->FS_values_and_derivatives(rhos, evdwl, dF_drho, ndensity);
-
- dF_drho_core = evdwl * dfcut + 1;
- for (DENSITY_TYPE p = 0; p < ndensity; ++p)
- dF_drho(p) *= fcut;
- evdwl_cut = evdwl * fcut + rho_core;
-
- // E0 shift
- evdwl_cut += basis_set->E0vals(mu_i);
-
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("dFrhos = ");
- for(DENSITY_TYPE p =0; pndensity;
- for (DENSITY_TYPE p = 0; p < ndensity; ++p) {
- //for rank=1 (r=0) only 1 ms-combination exists (ms_ind=0), so index of func.ctildes is 0..ndensity-1
- weights_rank1(func->mus[0], func->ns[0] - 1) += dF_drho(p) * func->ctildes[p];
- }
- }
-
- // rank>1
- func_ms_ind = 0;
- func_ms_t_ind = 0;// index for dB
- DOUBLE_TYPE theta = 0;
- for (func_ind = 0; func_ind < total_basis_size; ++func_ind) {
- ACECTildeBasisFunction *func = &basis[func_ind];
-// ndensity = func->ndensity;
- rank = func->rank;
- mus = func->mus;
- ns = func->ns;
- ls = func->ls;
- for (ms_ind = 0; ms_ind < func->num_ms_combs; ++ms_ind, ++func_ms_ind) {
- ms = &func->ms_combs[ms_ind * rank];
- theta = 0;
- for (DENSITY_TYPE p = 0; p < ndensity; ++p) {
- theta += dF_drho(p) * func->ctildes[ms_ind * ndensity + p];
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("(p=%d) theta += dF_drho[p] * func.ctildes[ms_ind * ndensity + p] = %f * %f = %f\n",p, dF_drho(p), func->ctildes[ms_ind * ndensity + p],dF_drho(p)*func->ctildes[ms_ind * ndensity + p]);
- printf("theta=%f\n",theta);
-#endif
- }
-
- theta *= 0.5; // 0.5 factor due to possible double counting ???
- for (t = 0; t < rank; ++t, ++func_ms_t_ind) {
- m_t = ms[t];
- factor = (m_t % 2 == 0 ? 1 : -1);
- dB = dB_flatten(func_ms_t_ind);
- weights(mus[t], ns[t] - 1, ls[t], m_t) += theta * dB; //Theta_array(func_ms_ind);
- // update -m_t (that could also be positive), because the basis is half_basis
- weights(mus[t], ns[t] - 1, ls[t], -m_t) +=
- theta * (dB).conjugated() * factor;// Theta_array(func_ms_ind);
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("dB(n,l,m)(%d,%d,%d) = (%f, %f)\n", ns[t], ls[t], m_t, (dB).real, (dB).img);
- printf("theta = %f\n",theta);
- printf("weights(n,l,m)(%d,%d,%d) += (%f, %f)\n", ns[t], ls[t], m_t, (theta * dB * 0.5).real,
- (theta * dB * 0.5).img);
- printf("weights(n,l,-m)(%d,%d,%d) += (%f, %f)\n", ns[t], ls[t], -m_t,
- ( theta * (dB).conjugated() * factor * 0.5).real,
- ( theta * (dB).conjugated() * factor * 0.5).img);
-#endif
- }
- }
- }
- energy_calc_timer.stop();
-
-// ==================== FORCES ====================
-#ifdef PRINT_MAIN_STEPS
- printf("\nFORCE CALCULATION\n");
- printf("loop over neighbours\n");
-#endif
-
- forces_calc_loop_timer.start();
-// loop over neighbour atoms for force calculations
- for (jj = 0; jj < jnum_actual; ++jj) {
- mu_j = elements[jj];
- r_hat = rhats[jj];
- inv_r_norm = inv_r_norms[jj];
-
- Array2DLM &Y_cache_jj = Y_cache(jj);
- Array2DLM &DY_cache_jj = DY_cache(jj);
-
-#ifdef PRINT_LOOPS_INDICES
- printf("\nneighbour atom #%d\n", jj);
- printf("rhat = (%f, %f, %f)\n", r_hat[0], r_hat[1], r_hat[2]);
-#endif
-
- forces_calc_neighbour_timer.start();
-
- f_ji[0] = f_ji[1] = f_ji[2] = 0;
-
-//for rank = 1
- for (n = 0; n < nradbasei; ++n) {
- if (weights_rank1(mu_j, n) == 0)
- continue;
- auto &DG = DG_cache(jj, n);
- DGR = DG * Y00;
- DGR *= weights_rank1(mu_j, n);
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("r=1: (n,l,m)=(%d, 0, 0)\n",n+1);
- printf("\tGR(n=%d, r=%f)=%f\n",n+1,r_norm, gr(n));
- printf("\tDGR(n=%d, r=%f)=%f\n",n+1,r_norm, dgr(n));
- printf("\tdF+=(%f, %f, %f)\n",DGR * r_hat[0], DGR * r_hat[1], DGR * r_hat[2]);
-#endif
- f_ji[0] += DGR * r_hat[0];
- f_ji[1] += DGR * r_hat[1];
- f_ji[2] += DGR * r_hat[2];
- }
-
-//for rank > 1
- for (n = 0; n < nradiali; n++) {
- for (l = 0; l <= lmaxi; l++) {
- R_over_r = R_cache(jj, n, l) * inv_r_norm;
- DR = DR_cache(jj, n, l);
-
- // for m>=0
- for (m = 0; m <= l; m++) {
- ACEComplex w = weights(mu_j, n, l, m);
- if (w == 0)
- continue;
- //counting for -m cases if m>0
- if (m > 0) w *= 2;
-
- DY = DY_cache_jj(l, m);
- Y_DR = Y_cache_jj(l, m) * DR;
-
- grad_phi_nlm.a[0] = Y_DR * r_hat[0] + DY.a[0] * R_over_r;
- grad_phi_nlm.a[1] = Y_DR * r_hat[1] + DY.a[1] * R_over_r;
- grad_phi_nlm.a[2] = Y_DR * r_hat[2] + DY.a[2] * R_over_r;
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("d_phi(n=%d, l=%d, m=%d) = ((%f,%f), (%f,%f), (%f,%f))\n",n+1,l,m,
- grad_phi_nlm.a[0].real, grad_phi_nlm.a[0].img,
- grad_phi_nlm.a[1].real, grad_phi_nlm.a[1].img,
- grad_phi_nlm.a[2].real, grad_phi_nlm.a[2].img);
-
- printf("weights(n,l,m)(%d,%d,%d) = (%f,%f)\n", n+1, l, m,w.real, w.img);
- //if (m>0) w*=2;
- printf("dF(n,l,m)(%d, %d, %d) += (%f, %f, %f)\n", n + 1, l, m,
- w.real_part_product(grad_phi_nlm.a[0]),
- w.real_part_product(grad_phi_nlm.a[1]),
- w.real_part_product(grad_phi_nlm.a[2])
- );
-#endif
-// real-part multiplication only
- f_ji[0] += w.real_part_product(grad_phi_nlm.a[0]);
- f_ji[1] += w.real_part_product(grad_phi_nlm.a[1]);
- f_ji[2] += w.real_part_product(grad_phi_nlm.a[2]);
- }
- }
- }
-
-
-#ifdef PRINT_INTERMEDIATE_VALUES
- printf("f_ji(jj=%d, i=%d)=(%f, %f, %f)\n", jj, i,
- f_ji[0], f_ji[1], f_ji[2]
- );
-#endif
-
- //hard-core repulsion
- DCR = DCR_cache(jj);
-#ifdef DEBUG_FORCES_CALCULATIONS
- printf("DCR = %f\n",DCR);
-#endif
- f_ji[0] += dF_drho_core * DCR * r_hat[0];
- f_ji[1] += dF_drho_core * DCR * r_hat[1];
- f_ji[2] += dF_drho_core * DCR * r_hat[2];
-#ifdef PRINT_INTERMEDIATE_VALUES
- printf("with core-repulsion\n");
- printf("f_ji(jj=%d, i=%d)=(%f, %f, %f)\n", jj, i,
- f_ji[0], f_ji[1], f_ji[2]
- );
- printf("neighbour_index_mapping[jj=%d]=%d\n",jj,neighbour_index_mapping[jj]);
-#endif
-
- neighbours_forces(neighbour_index_mapping[jj], 0) = f_ji[0];
- neighbours_forces(neighbour_index_mapping[jj], 1) = f_ji[1];
- neighbours_forces(neighbour_index_mapping[jj], 2) = f_ji[2];
-
- forces_calc_neighbour_timer.stop();
- }// end loop over neighbour atoms for forces
-
- forces_calc_loop_timer.stop();
-
- //now, energies and forces are ready
- //energies(i) = evdwl + rho_core;
- e_atom = evdwl_cut;
-
-#ifdef PRINT_INTERMEDIATE_VALUES
- printf("energies(i) = FS(...rho_p_accum...) = %f\n", evdwl);
-#endif
- per_atom_calc_timer.stop();
-}
\ No newline at end of file
diff --git a/lib/pace/ace_evaluator.h b/lib/pace/ace_evaluator.h
deleted file mode 100644
index 356ea52563..0000000000
--- a/lib/pace/ace_evaluator.h
+++ /dev/null
@@ -1,230 +0,0 @@
-/*
- * Performant implementation of atomic cluster expansion and interface to LAMMPS
- *
- * Copyright 2021 (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
- * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
- * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
- *
- * ^1: Ruhr-University Bochum, Bochum, Germany
- * ^2: University of Cambridge, Cambridge, United Kingdom
- * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
- * ^4: University of British Columbia, Vancouver, BC, Canada
- *
- *
- * See the LICENSE file.
- * This FILENAME 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 3 of the License, or
- * (at your option) any later version.
-
- * This program 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 this program. If not, see .
- */
-
-
-// Created by Yury Lysogorskiy on 31.01.20.
-
-#ifndef ACE_EVALUATOR_H
-#define ACE_EVALUATOR_H
-
-#include "ace_abstract_basis.h"
-#include "ace_arraynd.h"
-#include "ace_array2dlm.h"
-#include "ace_c_basis.h"
-#include "ace_complex.h"
-#include "ace_timing.h"
-#include "ace_types.h"
-
-/**
- * Basic evaluator class, that should accept the basis set and implement the "compute_atom" method using given basis set.
- */
-class ACEEvaluator {
-protected:
-
- Array2D A_rank1 = Array2D("A_rank1"); ///< 2D-array for storing A's for rank=1, shape: A(mu_j,n)
- Array4DLM A = Array4DLM("A"); ///< 4D array with (l,m) last indices for storing A's for rank>1: A(mu_j, n, l, m)
-
- Array1D rhos = Array1D("rhos"); ///< densities \f$ \rho^{(p)} \f$(ndensity), p = 0 .. ndensity-1
- Array1D dF_drho = Array1D("dF_drho"); ///< derivatives of cluster functional wrt. densities, index = 0 .. ndensity-1
-
- /**
- * Initialize internal arrays according to basis set sizes
- * @param basis_set
- */
- void init(ACEAbstractBasisSet *basis_set);
-
-public:
- // set of timers for code profiling
-
- ACETimer loop_over_neighbour_timer; ///< timer for loop over neighbours when constructing A's for single central atom
- ACETimer per_atom_calc_timer; ///< timer for single compute_atom call
-
-
- ACETimer forces_calc_loop_timer; ///< timer for forces calculations for single central atom
- ACETimer forces_calc_neighbour_timer; ///< timer for loop over neighbour atoms for force calculations
-
- ACETimer energy_calc_timer; ///< timer for energy calculation
- ACETimer total_time_calc_timer; ///< timer for total calculations of all atoms within given atomic environment system
-
- /**
- * Initialize all timers
- */
- void init_timers();
-
- /**
- * Mapping from external atoms types, i.e. LAMMPS, to internal SPECIES_TYPE, used in basis functions
- */
- Array1D element_type_mapping = Array1D("element_type_mapping");
-
-
- DOUBLE_TYPE e_atom = 0; ///< energy of current atom, including core-repulsion
-
- /**
- * temporary array for the pair forces between current atom_i and its neighbours atom_k
- * neighbours_forces(k,3), k = 0..num_of_neighbours(atom_i)-1
- */
- Array2D neighbours_forces = Array2D("neighbours_forces");
-
- ACEEvaluator() = default;
-
- virtual ~ACEEvaluator() = default;
-
- /**
- * The key method to compute energy and forces for atom 'i'.
- * Method will update the "e_atom" variable and "neighbours_forces(jj, alpha)" array
- *
- * @param i atom index
- * @param x atomic positions array of the real and ghost atoms, shape: [atom_ind][3]
- * @param type atomic types array of the real and ghost atoms, shape: [atom_ind]
- * @param jnum number of neighbours of atom_i
- * @param jlist array of neighbour indices, shape: [jnum]
- */
- virtual void compute_atom(int i, DOUBLE_TYPE **x, const SPECIES_TYPE *type, const int jnum, const int *jlist) = 0;
-
- /**
- * Resize all caches over neighbours atoms
- * @param max_jnum maximum number of neighbours
- */
- virtual void resize_neighbours_cache(int max_jnum) = 0;
-
-#ifdef EXTRA_C_PROJECTIONS
- /**
- * 2D array to store projections of basis function for rank = 1, shape: [func_ind][ndensity]
- */
- Array2D basis_projections_rank1 = Array2D("basis_projections_rank1");
-
- /**
- * 2D array to store projections of basis function for rank > 1, shape: [func_ind][ndensity]
- */
- Array2D basis_projections = Array2D("basis_projections");
-#endif
-};
-
-//TODO: split into separate file
-/**
- * Evaluator for C-tilde basis set, that should accept the basis set and implement the "compute_atom" method using given basis set.
- */
-class ACECTildeEvaluator : public ACEEvaluator {
-
- /**
- * Weights \f$ \omega_{i \mu n 0 0} \f$ for rank = 1, see Eq.(10) from implementation notes,
- * 'i' is fixed for the current atom, shape: [nelements][nradbase]
- */
- Array2D weights_rank1 = Array2D("weights_rank1");
-
- /**
- * Weights \f$ \omega_{i \mu n l m} \f$ for rank > 1, see Eq.(10) from implementation notes,
- * 'i' is fixed for the current atom, shape: [nelements][nradbase][l=0..lmax, m]
- */
- Array4DLM weights = Array4DLM("weights");
-
- /**
- * cache for gradients of \f$ g(r)\f$: grad_phi(jj,n)=A2DLM(l,m)
- * shape:[max_jnum][nradbase]
- */
- Array2D DG_cache = Array2D("DG_cache");
-
-
- /**
- * cache for \f$ R_{nl}(r)\f$
- * shape:[max_jnum][nradbase][0..lmax]
- */
- Array3D R_cache = Array3D("R_cache");
- /**
- * cache for derivatives of \f$ R_{nl}(r)\f$
- * shape:[max_jnum][nradbase][0..lmax]
- */
- Array3D DR_cache = Array3D("DR_cache");
- /**
- * cache for \f$ Y_{lm}(\hat{r})\f$
- * shape:[max_jnum][0..lmax][m]
- */
- Array3DLM Y_cache = Array3DLM("Y_cache");
- /**
- * cache for \f$ \nabla Y_{lm}(\hat{r})\f$
- * shape:[max_jnum][0..lmax][m]
- */
- Array3DLM DY_cache = Array3DLM("dY_dense_cache");
-
- /**
- * cache for derivatives of hard-core repulsion
- * shape:[max_jnum]
- */
- Array1D DCR_cache = Array1D("DCR_cache");
-
- /**
- * Partial derivatives \f$ dB_{i \mu n l m t}^{(r)} \f$ with sequential numbering over [func_ind][ms_ind][r],
- * shape:[func_ms_r_ind]
- */
- Array1D dB_flatten = Array1D