From 7740adb954686c1af7f7a7caea6bcc3339b2e35e Mon Sep 17 00:00:00 2001 From: wenyan4work Date: Sat, 17 Feb 2024 02:42:52 +0800 Subject: [PATCH] add torsionscan example --- .../advanced/geometric_torsionscan/.gitignore | 1 + .../geometric_torsionscan/geometric_scan.txt | 2 + .../geometric_torsionscan.py | 113 ++++++++++++++++++ .../geometric_torsionscan/plot_result.py | 74 ++++++++++++ 4 files changed, 190 insertions(+) create mode 100644 examples/advanced/geometric_torsionscan/.gitignore create mode 100644 examples/advanced/geometric_torsionscan/geometric_scan.txt create mode 100644 examples/advanced/geometric_torsionscan/geometric_torsionscan.py create mode 100644 examples/advanced/geometric_torsionscan/plot_result.py diff --git a/examples/advanced/geometric_torsionscan/.gitignore b/examples/advanced/geometric_torsionscan/.gitignore new file mode 100644 index 00000000..15b69c82 --- /dev/null +++ b/examples/advanced/geometric_torsionscan/.gitignore @@ -0,0 +1 @@ +**.png \ No newline at end of file diff --git a/examples/advanced/geometric_torsionscan/geometric_scan.txt b/examples/advanced/geometric_torsionscan/geometric_scan.txt new file mode 100644 index 00000000..62e90a24 --- /dev/null +++ b/examples/advanced/geometric_torsionscan/geometric_scan.txt @@ -0,0 +1,2 @@ +$scan +dihedral 21 20 5 3 -180.0 180.0 37 \ No newline at end of file diff --git a/examples/advanced/geometric_torsionscan/geometric_torsionscan.py b/examples/advanced/geometric_torsionscan/geometric_torsionscan.py new file mode 100644 index 00000000..e950e125 --- /dev/null +++ b/examples/advanced/geometric_torsionscan/geometric_torsionscan.py @@ -0,0 +1,113 @@ +# Copyright 2024 The GPU4PySCF Authors. All Rights Reserved. +# +# 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 . + +import time + +import pyscf +from pyscf import lib +from pyscf.geomopt.geometric_solver import optimize + +from gpu4pyscf.dft import rks + +lib.num_threads(16) + +# Enzalutamide +# https://www.echemi.com/cms/408806.html +# Enzalutamide is an androgen receptor (AR) inhibitor, +# which was launched in 2012 for the treatment of castrated prostate cancer and metastatic hormone-sensitive prostate cancer. +# coords in angstrom, generated by rdkit +atom = ''' +C -0.52489785 0.25103945 -2.93015020 +C -0.50680345 -0.64393976 -1.68199231 +C -1.91157268 -0.53966994 -1.07689036 +O -2.89596008 -1.08758214 -1.56376930 +N -1.87455596 0.28604347 0.02890000 +C -0.55589417 0.64601028 0.30659520 +S -0.03771290 1.85199711 1.35883674 +N 0.28457771 -0.08561063 -0.58428830 +C 1.69095398 -0.30104521 -0.44016185 +C 2.53850590 -0.42689013 -1.54537828 +C 3.90805921 -0.65833698 -1.38306826 +C 4.46679033 -0.77942388 -0.11616884 +C 3.63537497 -0.68825726 1.00120202 +C 2.26303281 -0.45702036 0.83390378 +C 5.91434343 -1.04523359 0.03306867 +O 6.48698930 -1.82113253 -0.71992937 +N 6.53927381 -0.36651782 1.06132882 +C 7.96812247 -0.45662800 1.23170297 +F 4.66727037 -0.72425525 -2.48600744 +C -3.05365343 0.62887069 0.75835844 +C -4.27570462 0.82269353 0.09079407 +C -5.46414648 1.14892461 0.78005884 +C -5.42793076 1.27078784 2.17688071 +C -4.22194909 1.06521787 2.85487200 +C -3.05316878 0.74497593 2.15329820 +C -6.60149929 1.59446376 2.93995702 +N -7.54355554 1.85428697 3.56412460 +C -6.74928175 1.36174974 0.01633060 +F -7.27014744 2.60177958 0.21630004 +F -6.60372258 1.23850678 -1.33208742 +F -7.71631245 0.47409741 0.36837460 +C -0.22537027 -2.11184985 -2.01127572 +H -0.73323150 1.29620611 -2.67177714 +H 0.42028392 0.23270869 -3.47581552 +H -1.30449790 -0.07030990 -3.63124037 +H 2.18898530 -0.32226984 -2.56341479 +H 4.03255894 -0.83066806 2.00353550 +H 1.64713582 -0.45176503 1.73036295 +H 6.06116959 0.41424928 1.49170789 +H 8.45082885 0.17225656 0.47895774 +H 8.22197650 -0.09849639 2.23196987 +H 8.30061119 -1.49102468 1.10983334 +H -4.31936045 0.74694915 -0.99665531 +H -4.18054201 1.14208036 3.94151726 +H -2.15934732 0.54471490 2.73896838 +H -0.25324848 -2.73066872 -1.10671674 +H -0.97426928 -2.51384795 -2.70351794 +H 0.75149209 -2.25342097 -2.48143483 +''' + +xc = 'HYB_MGGA_XC_WB97M_V' +bas = 'def2-tzvpp' +auxbasis = 'def2-tzvpp-jkfit' +scf_tol = 1e-10 +max_scf_cycles = 200 +screen_tol = 1e-14 +grids_level = 3 +mol = pyscf.M(atom=atom, basis=bas, max_memory=120000) + +mol.verbose = 1 +mf_GPU = rks.RKS(mol, xc=xc, disp=None).density_fit(auxbasis=auxbasis) +mf_GPU.grids.level = grids_level +mf_GPU.conv_tol = scf_tol +mf_GPU.max_cycle = max_scf_cycles +mf_GPU.screen_tol = screen_tol + +gradients = [] + + +def callback(envs): + gradients.append(envs['gradients']) + + +start_time = time.time() +mol_eq = optimize( + mf_GPU, + maxsteps=500000000, + constraints='geometric_scan.txt', # atom index is 1-based in this file + callback=callback) +print("Optimized coordinate:") +print(mol_eq.atom_coords()) +print(time.time() - start_time) diff --git a/examples/advanced/geometric_torsionscan/plot_result.py b/examples/advanced/geometric_torsionscan/plot_result.py new file mode 100644 index 00000000..61dacd8c --- /dev/null +++ b/examples/advanced/geometric_torsionscan/plot_result.py @@ -0,0 +1,74 @@ +import re + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.ticker import MaxNLocator + + +def parse_results(): + with open('scan-final.xyz', 'r') as file: + # Read all lines in the file + lines = file.readlines() + + # List to hold lines containing 'Dihedral' + dihedral_lines = [] + + # Loop through each line in the file + for line in lines: + # If 'Dihedral' is found in the line, add it to the list + if re.search('Dihedral', line, re.IGNORECASE): + dihedral_lines.append(line.rstrip()) + + records = [] + for dl in dihedral_lines: + cycle_txt, dihedral_txt, iteration_energy_txt = dl.split(';') + iteration_txt, energy_txt = iteration_energy_txt.split('Energy') + dihedral_degree = float(dihedral_txt.strip().split(' ')[-1]) + iteration_number = int(iteration_txt.strip().split(' ')[-1]) + energy_hartree = float(energy_txt.split(' ')[-1]) + cycle_number = int(cycle_txt.strip().split('/')[0].split(' ')[-1]) + records.append( + [cycle_number, dihedral_degree, iteration_number, energy_hartree]) + + return records + + +def plot_scan(records: list): + degrees = [x[1] for x in records] + iteration_numbers = [x[2] for x in records] + hartree = [x[3] for x in records] + + fig, ax1 = plt.subplots() + + color = 'tab:red' + ax1.set_xlabel('Dihedral (degree)') + ax1.set_ylabel('Energy (kcal/mol)', color=color) + ax1.plot(degrees, + 627.5096080305927 * (hartree - np.min(hartree)), + 'x-', + color=color) + ax1.tick_params(axis='y', labelcolor=color) + + ax2 = ax1.twinx() + + color = 'tab:blue' + ax2.set_ylabel('Geometric Optimization Iterations', + color=color) # we already handled the x-label with ax1 + ax2.plot(degrees, iteration_numbers, 'o-', color=color) + ax2.tick_params(axis='y', labelcolor=color) + ax2.yaxis.set_major_locator(MaxNLocator(integer=True)) + + plt.title('Energy and Geometric Optimization Iterations') + fig.tight_layout() # to adjust subplots to fit into the figure area. + + plt.savefig('scan_result.png', dpi=150) + + +def main(): + records = parse_results() + plot_scan(records) + for r in records: + print(r) + + +main()