shankar1729 / jdftx

JDFTx: software for joint density functional theory
http://jdftx.org
79 stars 49 forks source link

Coulomb Truncation Embed & within the margin of 5 bohrs from the truncation boundary error #333

Open AndrewTaehonKim opened 1 month ago

AndrewTaehonKim commented 1 month ago

Hello Shankar,

Thank you very much for the fantastic program you have built. I just have a quick question about a parameter in JDFTx and a potentially related error.

The question is regarding coulomb-truncation-embed.

coulomb-interaction Slab 001
coulomb-truncation-embed 0 0 7

I tried to look at the other forums you have commented on, the tutorials, and the documentation, but I am still not quite sure where to optimally place the coulomb-truncation-embed for my surface adsorption system.

Below are some places I have tried or am considering trying (see images below).

embbed

embbed2

Please let me know what the optimal place to put the coulomb-truncation-embed might be. I have been getting the following error frequently, and I am wondering if the poor computation speed and the following error may be due to that.

The output error is

Atom 74 lies within the margin of 5 bohrs from the truncation boundary.
Expand unit cell, or if absolutely sure, reduce coulomb-truncation-ion-margin.
End date and time: Sat Mar  2 21:07:03 2024  (Duration: 7-17:48:33.87)
Failed.

I thought the error was odd looking at the output xsf file because it appears that everything is well contained. Atom 74 is the second S atom in the ionpos file since there are a total of 78 atoms in the system. image

Here is the input, ionpos, and lattice data Input

include Na2S4@graphene_6x6x1-glyme.lattice
include Na2S4@graphene_6x6x1-glyme.ionpos
coords-type cartesian

coulomb-interaction Slab 001
coulomb-truncation-embed 0 0 5

ion-species GBRV/$ID_pbe.uspp
elec-cutoff 20 100

elec-smearing Fermi 0.01
kpoint-folding 1 1 1 

ionic-minimize \
    nIterations 200 \
    dirUpdateScheme FletcherReeves \
    energyDiffThreshold 1e-6 \
    knormThreshold 1e-4  #Threshold on RMS cartesian force

include ../SOLVENTS/glyme.in

dump-name Na2S4@graphene_6x6x1-glyme.$VAR
dump End State

Ionpos

# Ionic positions in cartesian coordinates:
ion C   28.01054201370563   0.002014969601156   0.13546520374716664 0
ion C   16.341968157882285  20.214705522223607  0.136370795694031   0
ion C   2.338532027351695   4.04241150539906    0.1303618477548163  0
ion C   23.34331318873597   0.002278716019722   0.14002246723956446 0
ion C   4.672620934470559   0.001554156700161   0.13205061600229584 0
ion C   11.674116326115334  20.215100669694028  0.11984749505502279 0
ion C   25.676530282256284  4.042904097616445   0.15356493008400918 0
ion C   18.6758513761297    16.172910511284943  0.16269558343134705 0
ion C   4.67177794833602    8.082036709586005   0.1125643862263992  0
ion C   18.676336142016087  0.002101964498535   0.1511803335222801  0
ion C   9.339850970925752   0.001738039847832   0.14151183696330705 0
ion C   9.341212652353349   16.17561329175046   0.09927930615555525 0
ion C   23.34294143972961   8.08492046601554    0.18324025457124904 0
ion C   21.009110272076708  12.128959849493356  0.1884533400045605  0
ion C   7.009085367604252   12.128528685817338  0.09726220467118907 0
ion C   14.00810547007721   0.001988816408382   0.1540088340121173  0
ion C   14.009120976896664  16.174984627607692  0.12993847789001123 0
ion C   7.005550759074101   4.041176716120078   0.11635089939380805 0
ion C   21.010101094155832  4.043302536041172   0.16867085070472854 0
ion C   16.344445762465487  12.130683361169906  0.16851880643879724 0
ion C   9.337643998902925   8.08172876398286    0.11493162272993374 0
ion C   16.34378301731489   4.042062713362387   0.16883618202480122 0
ion C   11.674085450205247  4.041278849785761   0.14219606057693568 0
ion C   11.675597583834264  12.13362224250672   0.13006825354940688 0
ion C   18.67807495615481   8.085220585829878   0.18576937438150765 0
ion C   7.007175797728255   20.216155157604103  0.11970932406816104 0
ion C   -6.996748017889804  20.214440105700564  0.15309468328938536 0
ion C   0.002354623270169   8.082966457809835   0.1574188223969486  0
ion C   14.012576459862723  8.083203697825988   0.1697128930694909  0
ion C   2.339067441062873   20.216465077180068  0.14622075717824856 0
ion C   -4.663670527593792  16.17249099218406   0.1785966420177587  0
ion C   2.336125814129785   12.127745765919109  0.1564175000031085  0
ion C   -2.332224039758668  12.1280302825893    0.18585986224579543 0
ion C   -2.32965935279241   20.215194179952768  0.15983032218419524 0
ion C   4.673933801454849   16.17774025221071   0.12535973450104976 0
ion C   0.002990588629148   16.174656174657997  0.16655442871845594 0
ion C   14.0082667201441    21.561374576259517  0.12977720083714672 0
ion C   2.338930816614256   1.348570645432834   0.1320097585278468  0
ion C   25.67710056297738   1.348988124751942   0.1419868952000769  0
ion C   16.342686998612663  17.520648868856785  0.1423020483360773  0
ion C   4.672307145938162   5.38821992196299    0.11605725655032373 0
ion C   21.01012349530289   1.349044430338556   0.1523948341468646  0
ion C   7.005937845436384   1.348033936218268   0.1280313892865088  0
ion C   11.674692189550257  17.521226022139913  0.1120577051666487  0
ion C   23.343784883173047  5.390651685198355   0.16961238986628047 0
ion C   18.67697254338436   13.478483594758318  0.1755959070335713  0
ion C   7.004963348139575   9.43211598652598    0.09753828326472913 0
ion C   16.34247207660575   1.348066366442404   0.15913017987028155 0
ion C   11.67352225634378   1.347802312891912   0.14654703450904805 0
ion C   9.340696631175405   13.481256291567632  0.09930387426130238 0
ion C   21.01048758662797   9.434228317035956   0.18954967082717644 0
ion C   9.340862584041272   21.56227531605612   0.12325120338741158 0
ion C   -9.330262385564788  21.561165527783945  0.14049142424325112 0
ion C   0.00447689367609    5.389148113018151   0.14593244579948106 0
ion C   14.011432130181314  13.480719160062115  0.14355688202397587 0
ion C   9.339051356617166   5.386992451363917   0.12107823949676799 0
ion C   18.677007069088155  5.390558998944344   0.17674105518834793 0
ion C   16.34572588460909   9.434637493571932   0.1791319995166809  0
ion C   11.674484029805171  9.432912355349794   0.1441944046974477  0
ion C   14.00942839582331   5.387478001177438   0.16273690443884803 0
ion C   4.673293824799766   21.561882602550448  0.13751922638322256 0
ion C   -6.99684147398532   17.520066852095106  0.1663466210473601  0
ion C   2.337134667876788   9.428913976861832   0.14338717810724688 0
ion C   -2.330217064543965  9.4326504527904 0.17856354836396804 0
ion C   -4.663349535264295  21.56148150508267   0.1523925737861731  0
ion C   7.008365318233412   17.523060541264897  0.1074506391357204  0
ion C   0.004487123343834   21.56167110349009   0.1541093181981985  0
ion C   -4.664147187494721  13.477588422800196  0.18804708273696846 0
ion C   4.67461980711176    13.479414560663225  0.12518806418918516 0
ion C   2.338448909138405   17.522455183042435  0.14563368911560914 0
ion C   -2.330712054942804  17.52079060661574   0.17005869677036856 0
ion C   0.000363098194401   13.478789474313889  0.17555922694800685 0
ion S   5.091970501871192   14.05098507736514   6.960415788839001   1
ion S   7.085763358365026   13.330850691324308  10.218518721987664  1
ion S   8.301157960236354   9.542429961599444   10.208687880127844  1
ion S   10.287593396578217  8.838620943050008   6.9425266949037585  1
ion Na  5.32495103010586    9.06963784764446    5.41995418971592    1
ion Na  10.060140479120443  13.838430842148847  5.445934273045367   1

Lattice

lattice Triclinic 28.0057986 28.0057986 32.1253 90 90 120

I have been getting this error for other systems as well, which seem well-contained in the lattice upon observing the output xsf file.

Please let me know what you think when you can.

Thank you very much for your time and guidance!

AndrewTaehonKim commented 1 month ago

Dear Shankar,

I apologize for the many questions as I am quite new to computational chemistry, especially ab initio QM using plane waves. It would be a huge help if you would provide some guidance on the coulomb-truncation-embed parameter and how you recommend to set it for a system of this type.

Thank you very much for your time and help.

shankar1729 commented 1 month ago

In general, put the embed center at the geometric center of the structure along truncated directions. The location along periodic directions is irrelevant, so in your diagrams, A and B are equivalent, and C and D are equivalent.

In that example, C or D are better because they are close to the geometric center, and that will put the trunctaion boundaries furthest from your topmost and bottom-most atoms.

Looking at your input, your topmost atom is ~ 10 bohrs and bottom-most is ~ 0 bohrs along z, so your 0 0 5 choice is optimal. When I run a dry run with your posted input, ionpos and lattice, I don't get an error. Can you confirm that you ran exactly these inputs when you got that atom out of boundary error?

I notice that you got the error after a week (based on the duration), likely after many ionic steps. Is your adsorbate drifting away during the run? If so, that could be the problem. Your setup is correct to start, but eventually your system is breaking your setup by becoming too large since the molecule becomes far from your surface. Visualize the trajectory of your molecule using the animated option of createXSF (it should work on a partial/failed output file too). (If you have used LAMMPS, this is similar to the lost atoms problem when using fixed boundaries.)

Also, looking at the atoms in your calculation, maybe you need dispersion corrections for your adsorbate to 'stick'? (This would need van-der-waals D3 in your input for example.

AndrewTaehonKim commented 1 month ago

Thank you for clarifying where to put the embed-center.

Using the createXSF, I was able to animate the geometries up to the error (animation file: https://drive.google.com/file/d/135q8EosWJZwtfBN_oBRmt5QV9YFBm9_4/view?usp=drive_link). It does appear that the adsorbate is floating upwards, drifting towards the boundary. In its final position in the animation, it still does not appear to have drifted towards the boundaries significantly (see image below of the last iteration before crash). image

Given the lattice c dimension is 32 bohrs, there should still be about around 16 bohrs left before it reaches the edge of the lattice. Would you be able to clarify why then I would still be getting the error: Atom 74 lies within the margin of 5 bohrs from the truncation boundary?

I will also look into implementing van-der-waals D3. I initially did not include it because it was expensive and some of my other sorbents like NiS2 should theoretically have very strong chemisorptive forces, but clearly for graphene, that might be necessary. Thank you for clarifying that.

shankar1729 commented 1 month ago

Okay, that was as expected then!

As for the reason for the error, the way the unit cell is rendered is misleading: remember that the box was periodic till you introduced the embed center at ~ 5 units. This breaks periodicity at half the box length above and below the center. Thus, the truncation boundary is roughly at 21 bohrs above the bottom in your cell, and your molecule only needs to get up to 16 bohrs to reach the margin of 5 bohrs.

Also, this is not an arbitrary choice with half the box being wasted. Remember that the wavefunctions and charge density of the atomic layer at z = 0 are extending several bohrs below, which means they are at the top end of the cell as well. That's why you can't allow the molecule to drift up all the way to the top of the box.

Finally, D3 disperson corrections should have negligible computational cost. It should at most add a few seconds per ionic step, and likely much less.