the fix of atomic position during AIMD

[luke419]

Hello,

First many thanks for your answers, Francois.
Can I fix the position of some atoms during AIMD simulation of qbox while other atoms are movable as usual?
if it is possible, would you let me know how I can do it?

Best regards,
Young

[fgygi]

Hello Young,
Yes, it can be done. You must use the constraint command documented in http://qboxcode.org/doc/html/usage/comm ... constraint .
To lock an atom in its current position, use a position constraint as e.g.

Code: Select all

 constraint define position cH1 H1

This will lock atom H1 in its current position during the simulation. Note that "cH1" is the name of the constraint, which can be used to identify the constraint e.g. if you want to delete the constraint. Note also that the constraint command must be used at the beginning of each MD simulation (it is not saved in the restart file).
If you want to constrain an atom to be at an arbitrary position, you can first use the move command to move the atom there.
Let me know if you have any questions.
Best,
Francois

[luke419]

Thank your for your answers, Francois.
Following your indication, I've run MD simulation of 128 water molecules in which 3 waters are fixed.
It worked well but after 164 cycles, it showed the following error message.

<electronic_charge ispin="0"> 1024.00000000 </electronic_charge>
<eigenvalue_sum> -343.38153289 </eigenvalue_sum>
<etotal_int> -2198.27336935 </etotal_int>
<electronic_charge ispin="0"> 1024.00000000 </electronic_charge>
<eigenvalue_sum> -343.38111858 </eigenvalue_sum>
<etotal_int> -2198.27378815 </etotal_int>
<electronic_charge ispin="0"> -nan </electronic_charge>
<eigenvalue_sum> -nan </eigenvalue_sum>
DoubleMatrix::potrf, info=1
DoubleMatrix::potrf, info=1
DoubleMatrix::potrf, info=1
DoubleMatrix::potrf, info=1

What would be wrong with it? Would you me me know how I can solve this?
I've attached the end part of the output file.

Best regards,
Young

[fgygi]

Hello Young,
Could you include the input file as an attachment?
Francois

[luke419]

Hello, Francois.

I've attached the input file.
Thank you for your helps.

Best regards,
Young

[fgygi]

Hello Young,
In the file md_128.i you need to use a separate command for each atom that is constrained to be fixed. The following input worked for 200 iterations:

Code: Select all

load gs_128.xml
constraint define position cO1 O1 
constraint define position cO2 O2 
constraint define position cO3 O3 
constraint define position cO4 O4 
constraint define position cH1 H1 
constraint define position cH2 H2 
constraint define position cH3 H3 
constraint define position cH4 H4 
constraint define position cH5 H5 
constraint define position cH6 H6 
constraint define position cH7 H7 
constraint define position cH8 H8 
set xc SCAN
set wf_dyn PSDA
set atoms_dyn MD
rseed
randomize_v 298
set thermostat BDP
set th_temp 298
reset_vcm
set ecut 65
#set dt 10
set dt 9.91735537190083
#set th_time 4958.67768595041
set iter_cmd compute_mlwf
set iter_cmd_period 10
run 200 5
save md_128_constraint_test.xml

The constraint command that you used has too many arguments, and should have been rejected with the following error message:

Code: Select all

[qbox] constraint define position con O1 O2 O3 O4 H1 H2 H3 H4 H5 H6 H7 H8
<cmd>constraint define position con O1 O2 O3 O4 H1 H2 H3 H4 H5 H6 H7 H8</cmd>
 Incorrect number of arguments for position constraint

which would cause the constraint command to be ignored.
I'm not sure if this could be the cause of the crash you observe. I used the above input file and could run 200 iterations successfully.

Note that you can put all the constraint commands in a file, e.g. constraints.i and then use the file name as a command in md_128.i.
Note also that water described with the SCAN functional requires a slightly elevated temperature, e.g. 330 K to have the properties of water at ambient conditions. See e.g. the paper J. Chem. Phys. 151, 164101 (2019); https://doi.org/10.1063/1.5124957

Also, the jump in <etotal> at the beginning of the MD run is likely due to the fact that you computed the ground state using PBE and then switched to SCAN for the MD run.

Best,
Francois