Development of reverse mapping system bridging dissipative particle dynamics and fragment molecular orbital calculation

We have developed a reverse mapping system to convert mesoscale structures generated by dissipative particle dynamics (DPD) simulations into nanoscale structures. This system is called DSRMS (DPD-based structure reverse mapping system) and is controlled by Python3 scripts using OCTA’s COGNAC program for DPD and molecular dynamics. The restored structures can be subjected to fragment molecular orbital calculations using the ABINIT-MP program for detailed nanoscale interaction analysis. Polyelectrolyte and lipid membranes have been used as illustrative example.

We have developed a reverse mapping system to convert mesoscale structures generated by dissipative particle dynamics (DPD) simulations into nanoscale structures.This system is called DSRMS (DPD-based structure reverse mapping system) and is controlled by Python3 scripts using OCTA's COGNAC program for DPD and molecular dynamics.The restored structures can be subjected to fragment molecular orbital calculations using the ABINIT-MP program for detailed nanoscale interaction analysis.Polyelectrolyte and lipid membranes have been used as illustrative example.© 2023 The Author(s).Published on behalf of The Japan Society of Applied Physics by IOP Publishing Ltd Supplementary material for this article is available online M olecular simulations have been used in the development of a wide range of materials science fields.A well-known difficulty with all-atom molecular dynamics (MD) simulations is the high dependence on the initial configuration due to polymer entanglement, etc when analyzing polymers or phase separation behavior.Therefore, coarse-grained (CG) simulations of mesoscale structures are expected to provide an alternative, efficient way to evaluate long-term behavior in particular.There are several CG models and methods tailored for specific applications, such as the Kremer-Grest model 1,2) to study the properties of polymer chains, dissipative particle dynamics (DPD 3) ) to predict phase separation behavior, and the Martini model, a general force field centered on biomolecules. 4,5)Conversely, CG methods cannot directly evaluate the interactions between atoms at the nanoscale.Socalled reverse mapping [6][7][8] has been developed to restore the nanoscale atomic structures. Rverse mapping methods have been improved [9][10][11][12][13][14][15][16][17][18][19][20][21] and recently machine learning techniques have been introduced.[22][23][24][25][26] The fragment molecular orbital (FMO) method was originated by Kitaura 27) and is now widely used in various fields such as biophysics and drug design.28,29) Large molecular systems consisting of more than tens of thousands of atoms are easily handled by FMO through parallel processing. Themost attractive feature of FMO calculation is the ability to perform detailed interaction analysis between fragments (or molecular units) in a non-empirical way from an application point of view. Takingdvantage of the reliability of FMO interaction calculations, we have established a method to calculate effective particle-particle interaction (or χ) parameters for DPD simulations 30) and systematized it as FCEWS (FMO-based Chi-parameter Evaluation Workflow System); 31) more recently, we have also introduced machine learning to increase the efficiency of parameter evaluation.32) The scheme of DPD with FMO-derived χ parameters was named FMO-DPD, and it has been successfully applied to several systems.[33][34][35][36][37] These applications may indicate that FMO-DPD with FCEWS is reliable as a multiscale method from the nanoscale to the mesoscale. Whe the reverse mapping procedure from the mesoscale to the nanoscale is developed in conjunction with the FMO calculations for interaction analysis, the bidirectionality of FMO-DPD is realized.Although the free energy differences are hardly calculated by DPD in contrast to MD, the set of FMO interaction energies 28,29) with restored structures could partially complement this issue.Accordingly, we have just developed DSRMS (DPD-based structure reverse mapping system) for this purpose.
][42] DPD simulations are then performed with the COGNAC module 43) of the J-OCTA suite of multiscale programs. 44,45)DSRMS is based on this framework and uses the following eight workflow steps controlled by Python3 scripts.
(i) Density conversion: to correlate the dimensionless quantities of DPD with real dimensions, a density conversion is performed.The conversion formula follows Groot's work. 46)In the case of a particle equivalent to four water molecules (120 Å 3 ), DPD length unit R c corresponds to 7.1 Å. (ii) Creation of an all-atom molecular information file corresponding to each particle: All-atom structures that are terminally capped are recorded, and atom IDs to be removed as terminals are also set (see Fig. S1).(iii) Collection of RESP charges: 47) the RESP charges of the registered all-atom molecules are obtained by using ABINIT-MP. 48,49)iv) Force field acquisition: the force field of the all-atom molecule of each registered segment is obtained.Currently Openbabel 50) is used and only the GAFF (General AMBER Force Field) set 51) can be acquired.(v) Atomization of DPD results: As shown in Fig. 1, the terminal atom of each segment molecule is rotated to align as closely as possible with the binding direction of the DPD particle.The segment molecule is then positioned and binds to the molecules assigned to the adjacent particle.At this stage, the mapped atomic structure is simply rotated and positioned without any relaxation.(vi) All-atom MD relaxation: the structure obtained in (iv) is relaxed using all-atom MD by COGNAC.The relaxation is performed until equilibrium is reached in the NPT ensemble.The progression from (v) to (vi) for a single molecule is shown in Fig. 2. (vii) Generation of fragmentation information: After reverse mapping, fragmentation information is generated for each type of molecule that appears in the system (see Fig. S2).(viii) Fragment Assignment and FMO Calculation: Fragments are assigned to the molecules cut to the size that can be calculated by FMO.FMO calculations and interaction analysis are finally performed with ABINIT-MP.
In conjunction with the reverse mapping by DSRMS, a couple of special settings are desired at the DPD stage.
(1) The molecular sizes contained in the particles are kept uniform.(2) The backbone potentials are imposed to ensure that the correct bond lengths and angles are maintained, and should ideally reflect the all-atom model of the target molecule.In this study, the bond length is defined using a harmonic-type potential between 1 and 2 particles, and the bond angle is defined by specifying the distance with a harmonic-type potential to 1-3 particles instead of an angle value. 35)1-2 particles are set to a force constant C = 160 kT and an equilibrium distance r e = 0.86 R c , and 1-3 particles are set to 80 kT and 1.72 R c .In addition, to maintain the sp 2 angle property (120 degrees), the r e of the particle containing the cis double bond was set to 1.49 R c .As an additional note, a constraint of approximately C = 500 kT is required to maintain a ring structure, such as cholesterol.
The reverse mapping process provides the following options for versatility.For substances such as water or hydrated ion, where multiple molecules exist within a particle, it is necessary to treat each as a separate molecule after reverse mapping.By specifying the "multi" mode for the molecular species, they can be recognized as multiple molecules based on bond lengths.For ions, presets including charge information can be invoked for water containing Na + and water containing Cl − .In addition, it is possible to map a given molecule to the center of multiple particles, rather than having one structure per particle.In the case of a multi-cyclic molecule (e.g.cholesterol), such a structure is broken down into multiple particles for DPD.However, when attempting to create an all-atom structure for each particle during reverse mapping, the bonds may become complex and the all-atom structure may collapse.Therefore, it is possible to solve this problem by restoring the structure to the center of multiple particles.
Figure 3 illustrates the application flow using FCEWS and DSRMS.Two examples of such a flow are the polymer electrolyte membrane (Nafion©) 34) and the lipid bilayer (POPC), 33) which were previously studied by FMO-DPD.Rather, DPD simulations were performed on a small-scale system of 1000 particles for Nafion and 5000 particles for POPC for demonstration purposes.For Nafion (see Fig. 4), an all-atom structure of 11 000 atoms was generated from the CG DPD structure (consisting of 1 000 particles), and the processed structure was subjected to an FMO-MP2 calculation using the 6-31 G(d′) basis set. 52)Water molecules (colored red) that directly coordinate and strongly interact with the sulfonic acid portion (yellow portion) at the interface between the electrolyte membrane polymer and water were identified.The averaged stabilization energy by IFIE was about −20 kcal mol −1 .For POPC (Fig. 5), an atomistic structure of about 60 000 atoms was obtained from the 5000 particle DPD result, and 6000 atoms near the water-lipid interface were extracted.The FMO calculation was then performed on this extracted structure in a similar manner.Water molecules (colored red) were identified that directly coordinate and strongly interact with the phosphate moiety Fig. 1.The process of mapping all-atom structures onto DPD particles.Each particle is assigned a pre-specified all-atom structure and, after being rotated to match the direction of the specified terminal atom (indicated by "T"), bonds are generated.In practice, the process is not performed on isolated molecules, but is applied collectively to the entire system.
110902-2 © 2023 The Author(s).Published on behalf of The Japan Society of Applied Physics by IOP Publishing Ltd (yellow portion) at the lipid membrane and water interface.The corresponding averaged interaction energy was approximately −27 kcal mol −1 .These evaluated interaction energies of both Nafion and POPC were reasonably close to the values calculated at the χ parameter determination stage using FCEWS, indicating mutual consistency between the nanoscale and mesoscale.
As an additional note, in these verification calculations the all-atom MD relaxation after reverse mapping was performed (with COGNAC 43) ) in a short time of 20 ps.A longer relaxation time may be required if there is a volume variation of the all-atom structures mapped onto CG particles.In addition, since multiple bulk and interface information can be obtained from a single point structure after MD calculations, where "IFIE" is an abbreviation for inter-fragment interaction energy. 29,38)The water (in red) forming clusters shows strong attractive interactions with the yellow part (sulfonic acid).110902-3 © 2023 The Author(s).Published on behalf of The Japan Society of Applied Physics by IOP Publishing Ltd it is desirable to analyze the interaction at multiple sites rather than at a specific site.Furthermore, by using multiple snapshots, more accurate interaction information can be obtained as a statistical average.
In this article, we have reported on the development of DSRMS, which allows mesoscale (or CG) DPD structures to be converted to nanoscale structures for interaction analysis by FMO calculations with ABINIT-MP. 38,39)The combination of FCEWS 30,31) and DSRMS can provide nano-meso bidirectionality for DPD simulations.The restored atomistic structures can be used as initial configurations for conventional all-atom MD, avoiding the problem of local minima.Besides FMO interaction analysis, some cutout sample structures are applicable to spectroscopic calculations with statistical averaging, such as X-ray absorption spectra of oxygen atoms in lipid membranes. 53)Finally, it should be noted that DSRMS will be released after additional usability improvements.

Fig. 2 .
Fig. 2.An illustration of reverse mapping and relaxation for a single molecule.The Roman numerals in the figure correspond to the scheme in the text.The Roman numerals (v) and (vi) in the figure correspond to the DSRMS workflow steps in the text.In practice, the process is not performed on isolated molecules, but is applied collectively to the entire system.

Fig. 3 .Fig. 4 .
Fig. 3. Illustrative scheme of FCEWS and DSRMS.Using FCEWS, nanoscale information is reflected in the χ parameters of DPD particles by FMO calculations.After reverse mapping by DSRMS, the recovered structure is subjected to FMO interaction analysis.(b) (c) (a)

Fig. 5 .
Fig. 5. Reverse mapping and FMO interaction analysis results for POPC.(a) DPD structure of 5000 particles.(b) Structure of 6000 atoms after reverse mapping of DPD and extraction near the water-lipid interface for FMO calculation.(c) Visualization of the interactions after FMO calculation.The color classification is the same as in Fig. 4. The water (in red) forming clusters shows strong attractive interactions with the yellow part (phosphoric acid).