Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2025 of Alexander Humeniuk.
3 :
4 : This file is part of the liquid_crystal plumed module.
5 :
6 : The liquid_crystal plumed module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The liquid_crystal plumed module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 : #include "core/ActionShortcut.h"
20 : #include "core/ActionRegister.h"
21 : #include "multicolvar/MultiColvarShortcuts.h"
22 :
23 : using namespace PLMD::multicolvar;
24 :
25 : namespace PLMD {
26 : namespace liquid_crystal {
27 :
28 : //+PLUMEDOC COLVAR NEMATIC_ORDER
29 : /*
30 : Calculate the nematic order parameter.
31 :
32 : The nematic order parameter S characterizes the orientational order of molecules
33 : and ranges from S=0 (isotropic) to S=1 (all molecular axes are perfectly parallel).
34 : Most liquids are isotropic, as there is no preferred direction, and have an order parameter
35 : close to 0. In liquid crystals, membranes and solids, molecules tend to align giving
36 : rise to order parameters closer to 1.
37 :
38 : $S$ is calculated from the distribution of the angles between the molecular axes ($\hat{u}_i$ for $i=1,\ldots,N$)
39 : and the nematic director $\hat{n}$,
40 : $$
41 : S = \frac{1}{N} \sum_{i=1}^N \left(\frac{3}{2} \cos^2(\theta_i) - \frac{1}{2} \right),
42 : $$
43 : with $\cos(\theta_i) = \hat{n} \cdot \hat{u}_i$.
44 :
45 : The nematic director depends on the distribution of the molecular axes
46 : and is computed as the eigenvector belonging to the largest eigenvalue
47 : of the $3 \times 3$ nematic order tensor,
48 : $$
49 : Q_{a,b} = \frac{1}{N} \sum_{i=1}^N \left(\frac{3}{2} u_{a,i} u_{b,i} - \frac{1}{2} \delta_{a,b} \right).
50 : $$
51 : Note that if only a single molecular axis is given (N=1), the nematic order parameter is always S=1.
52 :
53 : By adding a bias to the nematic order parameter, one can drive a liquid crystal from the
54 : isotropic to the nematic phase.
55 :
56 : The axis of a rod-like molecule is defined as the distance vector between two atoms,
57 : it points from the tail atom to the head atom.
58 :
59 : ```plumed
60 : # Assume there are three molecules with 20 atoms each.
61 : # In the first molecule the molecular axis vector points from atom 1 to atom 20,
62 : # in the second molecule it points from atom 21 to atom 40
63 : # and in the third from atom 41 to atom 60.
64 :
65 : # Compute nematic order parameter for the three molecules.
66 : S: NEMATIC_ORDER MOLECULE_STARTS=1,21,41 MOLECULE_ENDS=20,40,60
67 : PRINT FILE=colvar ARG=S
68 :
69 : # Add a bias to the nematic order parameter S.
70 : BIASVALUE ARG=S
71 : ```
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 : class NematicOrder : public ActionShortcut {
77 : public:
78 : static void registerKeywords(Keywords& keys);
79 : explicit NematicOrder(const ActionOptions&);
80 : };
81 :
82 : PLUMED_REGISTER_ACTION(NematicOrder,"NEMATIC_ORDER")
83 :
84 4 : void NematicOrder::registerKeywords(Keywords& keys) {
85 4 : ActionShortcut::registerKeywords( keys );
86 4 : keys.add("atoms","MOLECULE_STARTS","The atoms where the molecular axis starts.");
87 4 : keys.add("atoms","MOLECULE_ENDS","The atoms where the molecular axis ends.");
88 8 : keys.setValueDescription("scalar","the modulus of the average vector");
89 4 : keys.needsAction("CONSTANT");
90 4 : keys.needsAction("CUSTOM");
91 4 : keys.needsAction("DIAGONALIZE");
92 4 : keys.needsAction("DISTANCE");
93 4 : keys.needsAction("MATRIX_PRODUCT");
94 4 : keys.needsAction("TRANSPOSE");
95 4 : keys.needsAction("VSTACK");
96 4 : }
97 :
98 2 : NematicOrder:: NematicOrder(const ActionOptions& ao):
99 : Action(ao),
100 2 : ActionShortcut(ao) {
101 : // Fetch indices of atoms that define the tails and the heads of the molecular axes.
102 : std::vector<std::string> starts, ends;
103 2 : MultiColvarShortcuts::parseAtomList("MOLECULE_STARTS",starts,this);
104 4 : MultiColvarShortcuts::parseAtomList("MOLECULE_ENDS",ends,this);
105 :
106 2 : if( starts.size()!=ends.size() )
107 0 : error(
108 : "Mismatched numbers of atoms specified to MOLECULE_STARTS and MOLECULE_ENDS keywords. "
109 : "The molecular axes are specified by pairs of atoms."
110 : );
111 : // number of molecules with axes
112 2 : size_t number_of_axes = starts.size();
113 2 : std::string L = getShortcutLabel();
114 :
115 2 : if (number_of_axes == 1) {
116 : // Catch the edge case where there is only a single molecule. In this case the code further below
117 : // will not work, since VSTACK returns a (3,) vector instead of a (1,3) matrix.
118 :
119 : // If only a single molecular axis is given, the nematic order parameter is always S=1.
120 2 : readInputLine( L + ": CONSTANT VALUE=1.0");
121 : return;
122 : }
123 :
124 1 : std::string dlist = "";
125 5 : for(unsigned i=0; i<number_of_axes; ++i) {
126 : std::string num;
127 4 : Tools::convert( i+1, num );
128 8 : dlist += " ATOMS" + num + "=" + starts[i] + "," + ends[i];
129 : }
130 :
131 : // Calculate the lengths of the distance vectors
132 : // d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
133 2 : readInputLine( L + "_dvals: DISTANCE" + dlist );
134 : // Calculate the vectorial orientations of the molecules
135 : // dc: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
136 2 : readInputLine( L + "_dvecs: DISTANCE COMPONENTS " + dlist );
137 : // Convert the vectors into unit vectors
138 : // dux: CUSTOM ARG=dc.x,d FUNC=x/y PERIODIC=NO
139 : // duy: CUSTOM ARG=dc.y,d FUNC=x/y PERIODIC=NO
140 : // duz: CUSTOM ARG=dc.z,d FUNC=x/y PERIODIC=NO
141 2 : readInputLine( L + "_dux: CUSTOM ARG=" + L + "_dvecs.x," + L + "_dvals FUNC=x/y PERIODIC=NO");
142 2 : readInputLine( L + "_duy: CUSTOM ARG=" + L + "_dvecs.y," + L + "_dvals FUNC=x/y PERIODIC=NO");
143 2 : readInputLine( L + "_duz: CUSTOM ARG=" + L + "_dvecs.z," + L + "_dvals FUNC=x/y PERIODIC=NO");
144 : // Create a Nx3 matrix that contains all the unit vectors
145 : // u: VSTACK ARG=dux,duy,duz
146 2 : readInputLine( L + "_u: VSTACK ARG=" + L + "_dux," + L + "_duy," + L + "_duz");
147 : // uT: TRANSPOSE ARG=u
148 2 : readInputLine( L + "_uT: TRANSPOSE ARG=" + L + "_u");
149 : // Now compute the 3x3 matrix Q
150 : // Number of molecular axes.
151 : std::string N;
152 1 : Tools::convert( number_of_axes, N );
153 : // N: CONSTANT VALUE=N
154 2 : readInputLine( L + "_N: CONSTANT VALUE=" + N);
155 : // N x identity matrix
156 : // Id: CONSTANT VALUES=N,0,0,0,N,0,0,0,N NROWS=3 NCOLS=3
157 2 : readInputLine( L + "_Id: CONSTANT VALUES=" + N + ",0,0,0," + N + ",0,0,0," + N + " NROWS=3 NCOLS=3");
158 : // uTu: MATRIX_PRODUCT ARG=uT,u
159 2 : readInputLine( L + "_uTu: MATRIX_PRODUCT ARG=" + L + "_uT," + L + "_u");
160 : // Qsum: CUSTOM ARG=uTu,Id FUNC=((3*x-y)/2) PERIODIC=NO
161 2 : readInputLine( L + "_Qsum: CUSTOM ARG=" + L + "_uTu," + L + "_Id " + "FUNC=((3*x-y)/2) PERIODIC=NO");
162 : // divide by N
163 : // Q: CUSTOM ARG=Qsum,N FUNC=x/y PERIODIC=NO
164 2 : readInputLine( L + "_Q: CUSTOM ARG=" + L + "_Qsum," + L + "_N FUNC=x/y PERIODIC=NO");
165 : // Diagonalize Q
166 : // diag: DIAGONALIZE ARG=q
167 2 : readInputLine( L + "_diag: DIAGONALIZE ARG=" + L + "_Q");
168 : // Nematic order parameter is the largest eigenvalue of Q.
169 2 : readInputLine( L + ": CUSTOM ARG=" + L + "_diag.vals-1," + L + "_N FUNC=x PERIODIC=NO");
170 2 : }
171 :
172 : }
173 : }
|
