Report on Machine Learning Algorithm for Backmapping Macromolecules

Machine Learning Algorithm for Backmapping Macromolecules Project Title: Polymer Informatics Tools for Sustainable 3D Printing (PITS3D) Fellow: Dr. Petra Bačová Host Institution: Universidad de Cádiz, Spain Duration: 20 Months (January 2024– August 2025) 1. Executive Summary This report provides a brief summary of a deep learning-based methodology for reintroducing atomic detail into coarse-grained (CG) configurations of macromolecules. The approach employs a U-Net convolutional neural network (CNN) model and it has been successfully applied to poly(lactic acid) (PLA) systems [1], a key focus of the PITS3D project. More information about the methodology, performed simulations and developed tools can be found in ref. [1]. The CNN model is available on GitHub [2]. This part of the project was done in collaboration with The Cyprus Institute, more specifically with Eleftherios Christofi and Vagelis Harmandaris. 2. Machine Learning Algorithm Overview 2.1. Model Architecture The architecture of the Neural Network was inspired from the work done by Li et al. [3], where they introduced a general approach for backmapping CG macro-molecules by utilizing a “Pix2Pix" general-purpose conditional Generative Adversarial Network (cGAN) to perform an image-to-image translation. In this method the atomistic configurations were encoded from the XYZ vector components into Red-Green-Blue (RGB) values, treating backmapping as a super-resolution problem that maps low-resolution CG images to high-resolution atomistic ones. The method was employed, as an illustrative example, in CG models of homopolymers cis-1,4 polyisoprene melts. For simplicity, in our work [4,1] we utilize only the generator of the “Pix2Pix" cGAN, which is a U-Net based model. Also, instead of employing RGB values, the CNN was trained directly on atomic descriptors (namely, probability distribution functions of atomic bonds) that are capable of creating representative atomistic configurations in Cartesian coordinates of high molecular weight multi-component polymeric systems. The bonds among the connected pairs of atoms served as the target output of the artificial neural network, while the coordinates and the type of the CG particles served as input. After a number of experiments, where we examined the neural network’s behaviour for different depths and activation functions, we concluded that the neural network shown in Figure 1is reliable for this task. We utilize a U-net CNN based model, which consists of an encoder and a decoder network, with skip connections among them. For the encoder we stack five down-sample blocks which consist of a convolution layer with stride 2, a leaky ReLU activation function, and a batch-normalization layer. We note that we start with 64 filters and end up with 512. Then we pass the output of the encoder to the decoder network, where we have five up-sample blocks which consist of a transposed convolution layer with stride 2 and a ReLU activation function. For the first up-sample block we have a dropout layer with a rate of 0.5. We note that for the last layer of the network we have a transposed convolution layer with stride 1 and a tanh activation function, because we rescale the target output values in the interval [-1,1]. 1 Furthermore, for the training process we utilize mini-batch gradient descent with batches of size 64 and Adam optimization algorithm with an initial learning rate of 0.001, which was decreased down to 0.000001 by a factor of 8 once learning stagnates. The CNN was implemented in the open source TensorFlow 2 platform. [5] The computational time needed to train the model on a single NVIDIA Tesla V100-SXM2 GPU for 1000 epochs was around 24 hours. Figure 1: Schematic representation of the CNN used for the implementation of the method. Figure 5 from “Physics-Informed Deep Learning Approach for Reintroducing Atomic Detail in Coarse-Grained Configurations of Multiple Poly(lactic acid) Stereoisomers” by Eleftherios Christofi, Petra Bačová, and Vagelis A. Harmandaris. Licensed under CC BY 4.0. DOI: [10.1021/acs.jcim.3c01870] 2.2. Training Process and Loss Function To improve the accuracy of the networks we further employ prior knowledge both into the data preprocessing and into the generator objective. For this, the loss function of the neural network is augmented with several penalizing terms, based on prior knowledge of the physical properties of the molecular system under study, e.g., bond lengths, bond angles, and dihedral angles. We define the loss function as a linear combination of the loss function terms that penalize bond vectors (Lbv), bond lengths (Lbl), and dihedral angles Lda, via: ˜ L=λbvLbv +λblLbl +λbaLba +λdaLda (1) We treat the weights of the loss terms (i.e. values for λ= (λmae,λbl,λba,λda,λv0)) as hyperparameters. To find the optimum hyperparameter configuration of those λvalues we perform a quick grid search. To investigate structural similarity between the target and predicted atomistic configurations, we compute probability distribution functions of bond lengths, bond angles, and dihedral angles, as well as nonbonded pair radial distribution functions. For this, we define a metric parameter Ψ, as the total L1 difference of the probability densities between the target and prediction for all the distributions of bond 2 lengths, Pbl, bond angles, Pba, dihedral angles, Pda, and radial distribution function, g(r), through: Ψbl = nbl X i=1 |Pprediction bl,i (x)−Ptarget bl,i (x)|(2) Ψba = nba X i=1 |Pprediction ba,i (x)−Ptarget ba,i (x)|(3) Ψda = nda X i=1 |Pprediction da,i (x)−Ptarget da,i (x)|(4) In the above equations, nbl,nba, and nda denote the number of bond lengths, bond angles, and dihedral angles, respectively. Moreover, for the radial distribution function, Ψgr =|gprediction(r)−gtarget(r)|,(5) with rbeing the magnitude of the distance between non-bonded particles. Therefore, atoms directly participating to bond stretching, bond angle bending, and dihedral angle definitions are excluded from the radial distribution function. The sum of these quantities was denoted as Ψtotal. Lower Ψtotal values indicate higher prediction quality. 3. Application to Poly(lactic acid) (PLA) As an illustrative example of chiral biodegradable polymers, the algorithm was applied to amorphous PLA polymers, which may contain two types of stereoisomers: poly(L-lactide) acid (PLLA) and poly(D-lactide) acid (PDLA). The framework was tested on various model systems, from homopolymer stereoisomers of PLA to copolymers with randomly placed chiral centers (PDLLA). The objective was twofold: to develop an effective and versatile algorithm for backmapping chiral molecules at the all-atom scale and to design a tool capable of producing atomistic configurations of PLA with multiple molecular weights and compositions. The methodology uses an all-atom (AA) description. 3.1. Backmapping Procedure and Validation The procedure involves three main stages: preprocessing, training, and postprocessing. During preprocessing, the system’s CG coordinates and particle types serve as input, while probability distribution functions of atomistic bond vectors are collected as the target output. The training phase involves a U-Net CNN architecture described above. A crucial aspect is the loss function, which is augmented with physical penalizing terms for properties like bond vectors and bond lengths; preliminary tests showed that penalizing only these two terms (bond vectors and bond lengths) yielded the best results. Finally, postprocessing ensures the physical accuracy and desired stereochemistry of the generated atomistic configurations. The postprocessing described below includes a short sequence of atomistic simulations and customized codes to examine and correct the stereochemistry of the derived all-atom configurations. The codes are available open-access. [6] 1. Energy Minimization: After obtaining the backmapped structure, energy minimization is performed via a steepest descent algorithm. This is standard practice due to the ill-posed nature of the reverse problem. 2. Short Simulation: A very short simulation (0.1 ps) is run to slowly introduce excluded volume according to the selected force field. 3. Stereochemistry Checking and Correction: The trained model might misplace atoms, especially the hydrogen on the chiral atom, leading to incorrect stereochemistry. A custom-developed C code checks if the stereochemistry of each monomer corresponds to the desired sequence. If 3 not, the positions of the hydrogen atom and the methyl group are switched by a reflection matrix. This process mimics challenges in PLA synthesis where full stereochemistry control is difficult. 4. Additional Energy Minimization and Short Run: To address potential overlaps created by stereochemistry correction, an additional energy minimization step followed by another short run is performed. 5. Analysis Data Collection: Finally, data for analysis are collected from short MD simulations of a few nanoseconds (approximatelly 10 ns). The percentage of “wrong monomers” in the predicted configuration is generally low, allowing for effective correction through energy minimization. The backmapping model’s accuracy was rigorously validated by comparing its predicted atomistic structures against reference data derived from extensive atomistic molecular dynamics simulations. This validation focused on assessing both intramolecular and intermolecular structural deviations, primarily utilizing probability distribution functions of bond lengths, bond angles, and notably, dihedral angles. The results consistently showed that initial predictions had only minor deviations, which were largely eliminated after short molecular dynamics equilibration runs (point 5 above). Furthermore, the model effectively captured the atomic packing and local arrangements, demonstrated by strong agreement in intramolecular and intermolecular radial distribution functions and accurate reproduction of system densities. Overall, the validation process confirmed the model’s efficiency and efficacy in generating physically accurate atomistic configurations. This can significantly reduce computational demands for equilibration. References [1] E. Christofi, P. Bačová, and V. A. Harmandaris. Physics-Informed Deep Learning Approach for Reintroducing Atomic Detail in Coarse-Grained Configurations of Multiple Poly(lactic acid) Stereoisomers. Journal of Chemical Information and Modeling, 64(6):1853–1867, 2024. [2] E. Christofi. PLA Backmapping. https://github.com/SimEA-ERA/PLA-BackMap-CG. [3] W. Li, C. Burkhart, P. Polińska, V. A. Harmandaris, and M. Doxastakis. Backmapping CoarseGrained Macromolecules: An Efficient and Versatile Machine Learning Approach. The Journal of Chemical Physics, 153(4):041101, 2020. [4] E. Christofi, A. Chazirakis, C. Chrysostomou, M. A. Nicolaou, W. Li, M. Doxastakis, and V. A. Harmandaris. Deep Convolutional Neural Networks for Generating Atomistic Configurations of Multi-Component Macromolecules from Coarse-Grained Models. The Journal of Chemical Physics, 157(18), 2022. [5] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, et al. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, 2015. [6] P. Bačová. PLA Analysis Tools. https://github.com/pbacova/PLA_analysis_tools.git. 4