 So please open the Gromax practical 2 second part of practical in your browser. Here we will do already a real QMM simulations. But first I want to do a quick lecture recap of what you just saw. So how to deal with QMM, how to deal with point charges, how Gromax interface, Cp2k interface deals with them, how it deals with PBC and so on. Parallel binary conditions. So what is GIP? GIP is a QMM coupling scheme which is implemented in Cp2k. And how it is working. So basically this is your whole system now. You have QMPart and you have MMPart just green here. And you want somehow to include effect from the MMP point charges into the QMM system. How it is done in Cp2k? You expand each of your Gaussian electrostatic potential, which is just Columbia potential. It's 1 over R potential. You expand them in a sense like some of several Gaussians. And then you project that Gaussians onto the multigrid of your QM system like you do in the same way like you do in these QM basis functions. Gaussian basis functions. So this kind of unified model. And the only thing which remains is that our low term. So you have now several Gaussians. You expand your whole potential into Gaussians, some Gaussians plus some remaining part, but that remaining part is all smooth and small that it could be projected directly onto the courses possible grid. So this is how the GIP is working. So instead of Columbia potential you use now potential which consists of several Gaussian functions plus some remnants. And you project them onto the multigrid and you directly then use it from that multigrid into the QM. That multigrid is already embedded into your QM simulations. It's the same multigrid. Okay. Why it is important? Because, for example, now we want to do fully periodic QM. And how it's usually done? Not usually done, but how it's done in CP2K. It's also a very funny thing because it's a voice jeep and some additional scheme. So what you're doing? So you first project your MM, CP2K is doing well. It projects your MM point charges using the jeep onto the QM part. So like that. Now we have system like that in a small quantum box. So quantum box is usually much smaller than the whole system box, that green box. Then what you calculate in reality in CP2K, you calculate that thing. You calculate periodic QM box like that. And you converge using the quick step, the density of that box, the electron wave function, if you want to say it like that. But in reality, you calculate periodic system. But if you check, well, now, unless you have the same dimensions of the QM and the MM box are the same, you have a round radius here because your initial system was like that. So it was here, here and here boxes. But now it is looking like that, which is not correct. And to correct that in CP2K, CP2K uses so-called local scheme, which is from the paper, I guess. So what it is doing, it basically takes the density of the, takes the QM density, the multigrid, then it calculates interaction with a periodic cell with itself, basically, with periodic images of itself. Then it removes it completely, the interaction with the periodic cell. Then it shifts your QM system into the correct position in space. And then it recouples them back. So basically it calculates again the interaction between the QM system, periodic images of the QM system between themselves, but with the correct distance. So it do it like that. So it shift and then decouple, it's so-called decoupling and recapping scheme. So first you decouple your system image with its own periodic images. So you basically split the system and then you recouple them back with the correct periodicity, with the correct distance between them. You already have a correct QM system. Okay, it should be seen now, yeah? Okay, good. So this is how the QMM will work and this is how we reimplement some kind of engromax. So let's go to the tutorial. Let's make them... Let's do the real QMM system. The first system we will do will be the same student molecule, but now in a box of tips and few waters. Please do steps one to four of exercise three. And Rick, ask a question about the cell dimensions. So how is that generated? Yes, so basically there is two cells. I will show later. So the full cell is directly a M-M cell. It's directly inherited from the gromax. And QM cell dimensions is calculated automatically at the gromvp level, yes. So you don't need to worry about this. You still can adjust it. I will show again tomorrow how you can adjust it if you want. But it will create a large enough QM cell that your QM system will be inside that. Yes, in general this grom... Thomas asks if in general is gromax.cp2k I can use standard and dp files and add cp2k parameters flat or sound dp files and native batches. Yes, so in general just standard gromax.cp file just add this cp2k... gromax.cp2k interface stuff, which I already showed you, several parameters and it will automatically switch back. Make sure that for example your output that NST out options are set to a proper value or else you probably will never see any results because usually it's set up in classical simulations like some output structure every 10,000 steps or something like that, which will never happen. So it's always better to reduce these kind of things. But in general if you can run with the same setup or you have dissimulations, you can easily switch to QMM using just these parameters which I showed you in the cp2k... in that cp2k description. In case of pure cp2k QMM input files to the cell dimension, yes, in pure cp2k you need to set the cell dimensions manually and you need to consider that your cell... QM cell is large enough to contain all your QM atoms. The role of dispersion interaction between two atoms molecules in QMM, thanks in advance. Okay, so here we need to consider two things. So even three things. You have a dispersion interaction between the M-M and M-M atom. You have dispersion interaction between M-M and Q-M atom and you need dispersion interaction between Q-M and Q-M atom. So between M-M and M-M they are taken into account by the Wonderwall's parameters of Q-M and M-M for the field. Between Q-M and M-M atoms, they are also inherited from the Wonderwall's interaction inherited from the force field. So your Q-M atom still repulse M-M atoms also with Wonderwall's parameters. Between Q-M and Q-M atoms you need to use usually this dispersion correction. If you want to use such kind of things you need to use dispersion corrections. For example, Grimmy dispersion corrections. We shall begin discuss tomorrow and this is exactly what you want to do if you need them. But the role is there is a role. If you do not set up dispersion tracks, for example, between Q-M and M-M atoms, then you could see a collapse of your wave function at some points because your M-M point charge we just, well, there is, I think in the lectures I was discussed this kind of collapse problem when your electron from the Q-M part will be just tripped off of the Q-M part and just try to collapse onto the M-M point charge because the M-M point charge comes to close and there is no barrier to do that. And the dispersion interactions like Wonderwall they're adding this low, adding this repulsive behavior on the very low distances so that will never happen. Yeah, so there is also called polarization catastrophe some kind of points of polarity over hyper-polarization wave function but yeah, it is a known problem and whoever comes up usually Q-M atoms still repulse the M-M atoms with Wonderwall interactions from inherited from the force field. What are the disadvantages of setting up Q-M cell sites and M-M cell sites the same? Yeah, there is very one huge disadvantage that the bigger your Q-M cell is the more plane waves can fit inside that. The more plane waves can fit inside the Q-M box the harder your simulation will be. For example, if you increase your cell size Q-M cell size by a factor of two your simulations will be slower by a factor of 8 to 16 I guess. So you should always choose Q-M box size as large enough to contain all your Q-M atoms plus some additional space of course but it should not be big so that your simulation will be very slow. Yeah, that's why usually you do not set up them the same. You can in principle but that makes no sense at all because Q-M calculation you do not feasible. You will have too many plane waves. Yeah, basically. Okay, I think we can continue. So let's move on to our... Let's check how your Q-M input changes how CP2K input changes once you set up your system already with M-M point changes. So and first of all, we should go to that so basically do that less... Okay, please open with less that CP2K input file and go to that Q-M-M section in that and it should look like that. So and Q-M... In this section it defines the parameters to deal with Q-M-M, to deal with point charges inside the CP2K in general. So basically here is again the cell but this cell is not the same as in subsist cell and this is your Q-M cell. So Q-M cell is defined automatically as already have been asked it defined on the... at the ground level on the basis of dimensions of your Q-M system which will supply it to the grow box. So again it's fully pre-ordered. Here is eCouple, Gauss and use GIP-LIP 12 it is your Q-M-LIP coupling scheme and Gauss means using GIP, right? E-Coupling scheme and use GIP-LIP it's just a number of gaussians which you want to use in GIP expansion the more gaussians the harder will be or the harder will be your simulation but the more precise it will be. And regarding the question of Shudra Shan yes, the coupling is defined automatically it is some default behavior of the grow box CP2K interface but again you can change it tomorrow I can show you how you can change it. Also, Shiv asked could we use dispersion interaction energy to establish stability of complex so I'm importantly going to say again I would say no because it's too really too hard your complex, usually your complex is not defined by the dispersion interactions only it defined, I can say mostly but if your complex is not covalently bound complex then it defined mostly by the combative interactions so you need always to take into account all kind of interactions you have in the complex. Yes. Okay, let's continue with searching to the QMM so next what we see is that how we treat the periodicity in our QMM this is using of that level scheme so there is a periodic QMM so and there is some default parameters like multiple on and this is in force U2K that need to use level scheme to decouple and recouple you can also open CP2K out file with less CP2K output file with extension dot out the one that appears in your directory and you can see that in SCF cycles sometimes you will see this decoupling energy and recoupling energy this actually that level scheme so decoupling it means that interaction with this removing interaction with incorrect interaction with periodic box and recoupling is one shifted back okay and the next section you see for each kind of atoms indexes of the atoms which should be QMM so basically the list of indexes which should be QMM it's also automatically generated and for each kind of atoms which are seen in the QMM system from the other side if you check the MMM section you will see that everything here basically switched off so because Gromax fully handles description of the MMM region so you switch off the non-bonded forces you switch off evals summation you switch off PME so you don't do anything with the MMM you just inform that the CP2K that it doesn't need to do anything here and in the topology section there is a very important thing also automatically generated but here is the PDB file from which the CP2K picks up the point charges actually and it picks up from the extended beta field which starts at column 81 in the PDB file starting from that column starts the point charge in that file so it does not produce any connectivity and it just generates topologies later on because because interactions bonded interactions also will be treated by Gromax don't need to do that okay now yeah so if you again look into the steelben so opt input you will see that only the atoms which belongs to steelben is marked as the QM atoms everything else will be treated as the MMM atoms and there is how the PDB file that PDB file looks like generated so basically here starts here is a splitting between two residue names QM and MMM I think it's self-explanatory so the QM is QM atoms and MMM atoms and here starts the point charges so you see that all QM atoms have a zero point charge it's verified and all MMM atoms have some point charge here for example you can see that these three atoms forms a deep city water pretty remarkable point charges of deep city water okay but anyway let's see what is the outcome of the exercise 3 I think it should finish and you can do this remaining part of the exercise 3 which is from 5 to ah no it's not remaining part it's from 5 to 7 that energy minimization and if you do that if you again extract the potential energy you will see something like that and this side just rendered your trajectory of that optimization and you see that your energy so there is much more steps that you have done but yeah you'll see how the energy minimizes ah again so yeah I guess now you can do the step 7 to 9 so now you instead of the optimization you again run the molecular dynamics using the NVT ensemble actually NPT is also working so in case you wondering can use a pressure coupling yes you can use pressure coupling ah yeah but for that so just use this one is so Thomas asking is 1 femtosecond type such standard focal length calculations wow it's a very good question what time step you should use normally I would say yes in most cases 1 femtosecond type step is kind of good for the QMM molecular dynamics except of several important things so you should know how you choose a time step for example in normal classical micro dynamics you should choose your time step according to the vibrations of your hydrogen atoms because it means that vibrations of your bonds some of your bonds should be period of that period of that vibration should be large enough that your molecular dynamics is sampled correctly and or you can use constraints instead of the bonds yeah so there is a constraints for example in Gromach since you set up constraints you can use even higher time steps like not 1 femtosecond but 2 to 3 femtoseconds but for that you need to use constraints because so else your hydrogen your hydrogen containing bonds will be vibrated too fast with respect to the time step but in most cases for QMM it is kind of okay to have 1 femtosecond time step if you see if you want to study the problem connecting to the which connects somehow to the vibrations of hydrogen for example if you want to transfer to study transfer of hydrogen from one group to another of hopping from one water to another or between the protein inside the protein between several amino acids then I would suggest you to reduce the time step to 0.5 femtosecond to a little lower Merth is asking in other bioexcel can we workshop the first part where you can define the the anthropology is that necessary or are you conducting say yeah so idea is that here it is all done automatically I guess yeah so the interface stays care of all that kind of parameters so you don't need to do anything by your hand in adjusting parameters but I see what you are implying that usually the hydrogen does not have any repulsion parameters so it hasn't any repulsion the hydrogen does not have any wonderful parameters yes it's true and it's true that it could be a problem I would suggest you to always check if your hydrogen atom from your water or from another part tries to collapse onto your part if it is the case then yes better to add to that hydrogen repulsive energy potential because it's an artificial thing yes for that yes Thomas Merth actually additionally we are also not providing any topology files to CP2K also this interface so you see the charge of atoms through the PDP file yes you are correct so that's the advantage of the interface that you don't need to provide almost an apology to the CP2K the Gromax is doing it for you and the Gromax takes care of the MM part of the MM topology the CP2K only see the QM atoms and point charges around that QM atoms that's all he see and everything else is treated by Gromax Umesh asks how frequently does one dump the coordinates into the output it's up to you clearly I can say just be warned that well practically usually I am doing this each step in QMM setup because you cannot have as many steps as you are doing this usually but it's up to you how many and how you want to deal with the gigabytes or even bigger files yeah I would say I would say as frequently as you can afford yes which makes sense yes it just affects this output does not affect like at all the speed of your QMM simulation because anywhere the QM will be a limiting factor bottleneck of your simulation the output is not here the limiting factor like a normal classical so you can do as frequently as you want as much as you can afford or analyze afterwards okay I suppose most of the people done is at parts 7 to 9 and here I show what should be as a result so if you you can also download trajectory and render it if you want and here is what you see in the temperature actually that's what you usually see what happens when you switch your system from MMM to QMM so basically here initial frame was a frame taken from the classical and trajectory and that's what happens usually when you instantaneously switch from MMM to QMM so you see a huge jump in temperature because your quantum system quantum part starts relaxing from the force field bond's length angles and so on and you see a huge jump in temperature because your energy starts dropping and you see increase of kinetic energy your potential drops and your kinetic energy increases and you see that temperature jump and then the thermostat starts to equilibrate your system so that's why I told you at some point before that you always need to after the trajectory you always need to relax your system for some time with QMM and then with the ensemble just to make sure that your system is well-equilibrated for QMM I guess so that was the idea of that that's the idea of that exercise so to see that what's happening when you switch your topology from QMM to QMM instantaneously on one step okay so let's move to the next exercise and in that exercise I want to show you what will look like if you set up really huge protein system and for that system I've used one of the proteins which I work in this so it's called phytochrome protein it's really huge protein it's for active protein it's contained on standard residue chromophore here in the pockets there is two of them here and there and the objective for the original work was to using the umbrella sampling to make the isomerization profile of that chromophore from that position to that position so it's basically a ring flip from that to that and that here will be a chromophore plus several amino acid residues around so it will be really big system it's round 85 quantum atoms something like that and for the QMM method we'll still use APB functional and amount 4 will be the ampere 0,3 so yeah do the steps 1,2,3 you can also download and render your the system in any software package you want but it's already really a big protein and you will notice how slower the simulation will be when you run it also it contains I think around 70 to 80 thousand mm point charges if I'm not mistaken but it's kind of a big system so there is several questions here so Merth asking what happens to total mm charts afterwards do we need to modify the charts in order to remain the same yes so that's one of the problems which I'm now trying to solve in addition to the interface so if you look into the output of the ground pp it notifies you at some point that your total QMM plus mm charts should be in principle the same as it was before in mm system however especially in a protein case usually your QMM plus mm charts will be not zero and even not the values that you want it to be so what to do with that for now the interface warns you that you have access charts and you need manually to spread it along the... you don't need usually it's not so important it's usually a small charge but you need to spread it along the nearby mm atoms to to have the same mm charts as it was before I'm still thinking how it can be automatized because it's not so easy problem but if your charge is small enough so if it's not like minus 1 or minus 2 X charge then you could simply use non zero charge much yes that from my experience I usually do not do that unification it works quite fine so if you distribute... what if you distribute the remaining circle of charge the same residue proportionally on the negative sense but yes so that's what you need to do and that's what if you look into the ground pp output that's what interface suggest you to do so you just need to... usually how it is done it is done that you spread that excess charge onto the mm atoms nearby to the key part in order to keep total charge the same as it was before in mm system so Shudrashan asks I suppose simulation can be restart the same way as it do classical dissimulation gromax so you can use restart cpt files generated by the gromax you can restart from the cpt checkpoint files or you can regenerate and restart or you can use cpt files however I would suggest you to then put additional option to the empty run I think it should be minus cpt minus cpt so the time interval between the gromax writes that cpt files it's in real time so they put it to something like minus cpt 0.1 or something like that so it will write cpt file on each step of your KMM simulation that's what I suggest usually to do but yes you can use cpt restart so in addition I can divide surplus charge among the residues remaining in the mm and subtract from the minuda I would suggest to divide this charge around the onto the mm atoms which are really nearby to the QM system to the broken QM bond so you have a bond which are broken between QM and MMParts that's why I suggest to put them on to the first atom on to the next atoms and you can define them evenly is a standard practice to perform might you ask is it possible I know you mentioned that MMPT is possible but you didn't recommend so just curiously why so yes typically the NBT ensemble I would prefer to do the NBT ensemble because NBT is just so I mean it makes sense because in practice all ensembles should give you the same results in a long distance trajectories like this ergoidal hypothesis tells you about that but with pressure coupling what could happen that you start changing your because how usually the pressure coupling is done in all microdynamic software you start changing the dimensions of your box and of course it's not feasible to change the dimensions of your you can change freely the dimensions of your MMPox pretty easily it is very well known algorithms but you cannot really change the dimensions of your QMbox effectively because once you start changing the dimensions of your QMbox you start changing your plane wave basis and you cannot do that because you start changing your QM system how it is implemented for example in Gromax in CP2K QMM dynamics in normal they always use the same QMbox for calculating QMPox but they also have a fake QMbox to pressure coupling so it makes your life more miserable so I've also made that easier so now the pressure coupling is fully implemented in Gromax so Gromax take care of that and how it takes care it scales only the MMPox it does not touch the QMbox at all so in a sense you don't have a pressure coupling of your QM part but except of that the NPT works quite well you can also try that yourself if you want yeah but I would suggest for QMM especially in CP2K to just stick to the NPT ensemble you can do preliminary NPT equilibration with normal MADECOR dynamics and then switch to NPT and use NPT yes, okay so I think most of the people done with parts steps one to three and let's continue you can open with less the phytogram input file and what you see will be in QMM appears additional section called link several link sections and in that link section what it is it is basically it defines the broken bonds between the QMM and MMPart like in that case of protein for example here is our QMPart and what you see here is that here, here, here and here your carbons doesn't have full balances and of course this is not allowed in the quantum chemistry dangling bonds kind of this will be a very bad for the QM so what in reality is going on you're setting up that you're saying between these two atoms we have a broken covalent bond and it informs the CP2K that it needs to generate here additional fake atom which is not which is not used in any real simulations but it just caps additional hydrogen atom will cap this free bond and after the QM calculations will be removed back so the force from that QM atom will just go on to the nearby QMM atom will be redistributed and this atom will be removed so but it just informs CP2K that it need to cap this broken bonds by fake atom it's called link atom I guess it was also in the lecture okay so in a way the result of this simulation I could show you so basically we cannot make it further because it's really hard simulation it takes months to simulate fully but in reality what you can get on the output is something like that so it's again the umbrella sampling of that reaction with that reaction coordinate so it's very smooth kind of profile but it took almost 40 picoseconds so there was like 20 again frames each several degrees you have a frame and it is simulated each frame for 30 to 40 picoseconds so it's really long simulation so one simulation took approximately couple of weeks of real computing time so it's really hard problem but that's what really you are usually dealing with in real life and that's why also important very important things that's why you need a supercomputers you cannot do this single kind of simulation on your laptop you need supercomputers and that's why the Arno me and Holly also trying to make you getting to know how to use them okay I think this is end for today practical so you can continue your continue doing whatever you did not done so I think now it still should work so thanks for your attention we are ready for some questions so I still read the question from Victoria does code run mostly on CPU and GPU this code is now run only on the CPU we do not use GPU however we have a plans to also make it work with GPU in the future I guess yes Arno according to what we have certainly I mean so Arch2o I certainly right now we are not using it with GPUs but in principle should be possible we know Gromax of course takes very good advantage of GPUs, CP2K can also, the actual advantage that you will get will depend on the system but that is something that we will providing guidance on the future so tomorrow Holly will talk a bit about performance and it will be mostly about performance actually I think almost totally on CPU the idea is that there is a best practice guide which you will refer to and we will in future expand that to include guidance on how to get performance on CPU plus GPU can't we still utilize Gromax to use GPU yeah you can, I mean it will be totally fine but there is no point to do that because your bottleneck in the QMM calculations always the QMPart regardless so regardless of anything the Gromax I mean the MMM calculation and integration step takes less than 0.01% of the total time so why it should use GPU is not clear for me the good thing will be once we can make the interface once we could make in the interface use of CPU by CP2K if CP2K can efficiently use GPU it will be much more increasing in speed than for Gromax, yeah you can so I mean the answer to the math questions yes you can but there is not big point to do that for now but we yeah not next to nothing you can even see a degradation of the performance because the time to submit of exchange between a CPU and GPU will be more than the actual calculation time on the GPU it will be not a huge degradation again around 0.0% but there is no point for now to do that once we can get into the GPU using GPU for CP2K simulations it will be much more clear why you should do that okay yeah so thanks very much Dimitri good sessions