-
Notifications
You must be signed in to change notification settings - Fork 13
/
module_QMMM.py
1424 lines (1233 loc) · 72.8 KB
/
module_QMMM.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import copy
import time
import numpy as np
import math
import shutil
#functions related to QM/MM
import ash.modules.module_coords
from ash.modules.module_coords import Fragment, write_pdbfile
from ash.functions.functions_general import ashexit, BC,blankline,listdiff,print_time_rel,printdebug,print_line_with_mainheader,writelisttofile,print_if_level
import ash.settings_ash
#QM/MM theory object.
#Required at init: qm_theory and qmatoms. Fragment not. Can come later
#TODO NOTE: If we add init arguments, remember to update Numfreq QMMM option as it depends on the keywords
class QMMMTheory:
def __init__(self, qm_theory=None, qmatoms=None, fragment=None, mm_theory=None, charges=None,
embedding="Elstat", printlevel=2, numcores=1, actatoms=None, frozenatoms=None, excludeboundaryatomlist=None,
unusualboundary=False, openmm_externalforce=False, TruncatedPC=False, TruncPCRadius=55, TruncatedPC_recalc_iter=50,
qm_charge=None, qm_mult=None, dipole_correction=True):
module_init_time=time.time()
timeA=time.time()
print_line_with_mainheader("QM/MM Theory")
#print(BC.WARNING,BC.BOLD,"------------Defining QM/MM object-------------", BC.END)
#Check for necessary keywords
if qm_theory is None or qmatoms is None:
print("Error: QMMMTheory requires defining: qm_theory, qmatoms, fragment")
ashexit()
#Defining charge/mult of QM-region
self.qm_charge=qm_charge
self.qm_mult=qm_mult
#Indicate that this is a hybrid QM/MM type theory
self.theorytype="QM/MM"
#External force energy. ALways zero except when using openmm_externalforce
self.extforce_energy=0.0
#Linkatoms False by default. Later checked.
self.linkatoms=False
#Whether we are using OpenMM custom external forces or not
#NOTE: affects runmode
self.openmm_externalforce=openmm_externalforce
#Theory level definitions
self.printlevel=printlevel
self.qm_theory=qm_theory
self.qm_theory_name = self.qm_theory.__class__.__name__
self.mm_theory=mm_theory
self.mm_theory_name = self.mm_theory.__class__.__name__
if self.mm_theory_name == "str":
self.mm_theory_name="None"
print("QM-theory:", self.qm_theory_name)
print("MM-theory:", self.mm_theory_name)
#If fragment object has been defined
#This probably needs to be always true
if fragment is not None:
self.fragment=fragment
self.coords=fragment.coords
self.elems=fragment.elems
self.connectivity=fragment.connectivity
# Region definitions
self.allatoms=list(range(0,len(self.elems)))
print("All atoms in fragment:", len(self.allatoms))
#Sorting qmatoms list
self.qmatoms = sorted(qmatoms)
if len(self.qmatoms) == 0:
print("Error: List of qmatoms provided is empty. This is not allowed.")
ashexit()
self.mmatoms=listdiff(self.allatoms,self.qmatoms)
# FROZEN AND ACTIVE ATOM REGIONS for NonbondedTheory
if self.mm_theory_name == "NonBondedTheory":
#NOTE: To be looked at. actatoms and frozenatoms have no meaning in OpenMMTHeory. NonbondedTheory, however.
if actatoms is None and frozenatoms is None:
#print("Actatoms/frozenatoms list not passed to QM/MM object. Will do all frozen interactions in MM (expensive).")
#print("All {} atoms active, no atoms frozen in QM/MM definition (may not be frozen in optimizer)".format(len(self.allatoms)))
self.actatoms=self.allatoms
self.frozenatoms=[]
elif actatoms is not None and frozenatoms is None:
print("Actatoms list passed to QM/MM object. Will skip all frozen interactions in MM.")
#Sorting actatoms list
self.actatoms = sorted(actatoms)
self.frozenatoms = listdiff(self.allatoms, self.actatoms)
print("{} active atoms, {} frozen atoms".format(len(self.actatoms), len(self.frozenatoms)))
elif frozenatoms is not None and actatoms is None:
print("Frozenatoms list passed to QM/MM object. Will skip all frozen interactions in MM.")
self.frozenatoms = sorted(frozenatoms)
self.actatoms = listdiff(self.allatoms, self.frozenatoms)
print("{} active atoms, {} frozen atoms".format(len(self.actatoms), len(self.frozenatoms)))
else:
print("active_atoms and frozen_atoms can not be both defined")
ashexit()
#print("List of all atoms:", self.allatoms)
print("QM region ({} atoms): {}".format(len(self.qmatoms),self.qmatoms))
print("MM region ({} atoms)".format(len(self.mmatoms)))
#print_time_rel(timeA, modulename="Region setup")
timeA=time.time()
#print("MM region", self.mmatoms)
blankline()
#TO delete
#List of QM and MM labels
#self.hybridatomlabels=[]
#for i in self.allatoms:
# if i in self.qmatoms:
# self.hybridatomlabels.append('QM')
# elif i in self.mmatoms:
# self.hybridatomlabels.append('MM')
else:
print("fragment= keyword has not been defined for QM/MM. Exiting")
ashexit()
#Setting QM/MM qmatoms in QMtheory also (used for Spin-flipping currently)
self.qm_theory.qmatoms=self.qmatoms
#Setting numcores of object.
#This will be when calling QMtheory and probably MMtheory
#numcores-setting in QMMMTheory takes precedent
if numcores != 1:
self.numcores=numcores
#If QMtheory numcores was set (and QMMMTHeory not)
elif self.qm_theory.numcores != 1:
self.numcores=self.qm_theory.numcores
#Default 1 proc
else:
self.numcores=1
print("QM/MM object selected to use {} cores".format(self.numcores))
#Embedding type: mechanical, electrostatic etc.
self.embedding=embedding
print("Embedding:", self.embedding)
#Whether to do dipole correction or not
#Note: For regular electrostatic embedding this should be True
# Turn off for charge-shifting
self.dipole_correction=dipole_correction
#if atomcharges are not passed to QMMMTheory object, get them from MMtheory (that should have been defined then)
if charges is None:
print("No atomcharges list passed to QMMMTheory object")
self.charges=[]
if self.mm_theory_name == "OpenMMTheory":
print("Getting system charges from OpenMM object")
#Todo: Call getatomcharges directly or should that have been called from within openmm object at init ?
#self.charges = mm_theory.getatomcharges()
self.charges = mm_theory.charges
elif self.mm_theory_name == "NonBondedTheory":
print("Getting system charges from NonBondedTheory object")
#Todo: normalize charges vs atom_charges
self.charges=mm_theory.atom_charges
else:
print("Unrecognized MM theory for QMMMTheory")
ashexit()
else:
print("Reading in charges")
if len(charges) != len(fragment.atomlist):
print(BC.FAIL,"Number of charges not matching number of fragment atoms. Exiting.",BC.END)
ashexit()
self.charges=charges
if len(self.charges) == 0:
print("No charges present in QM/MM object. Exiting...")
ashexit()
#Flag to check whether QMCharges have been zeroed in self.charges_qmregionzeroed list
self.QMChargesZeroed=False
#CHARGES DEFINED FOR OBJECT:
#Self.charges are original charges that are defined above (on input, from OpenMM or from NonBondedTheory)
#self.charges_qmregionzeroed is self.charges but with 0-value for QM-atoms
#self.pointcharges are pointcharges that the QM-code will see (dipole-charges, no zero-valued charges etc)
#Length of self.charges: system size
#Length of self.charges_qmregionzeroed: system size
#Length of self.pointcharges: unknown. does not contain zero-valued charges (e.g. QM-atoms etc.), contains dipole-charges
#self.charges_qmregionzeroed will have QM-charges zeroed (but not removed)
self.charges_qmregionzeroed = []
#Self.pointcharges are pointcharges that the QM-program will see (but not the MM program)
# They have QM-atoms zeroed, zero-charges removed, dipole-charges added etc.
#Defined later
self.pointcharges = []
#Truncated PC-region option
self.TruncatedPC=TruncatedPC
self.TruncPCRadius=TruncPCRadius
self.TruncatedPCcalls=0
self.TruncatedPC_recalc_flag=False
self.TruncatedPC_recalc_iter=TruncatedPC_recalc_iter
if self.TruncatedPC is True:
print("Truncated PC approximation in QM/MM is active.")
print("TruncPCRadius:", self.TruncPCRadius)
print("TruncPC Recalculation iteration:", self.TruncatedPC_recalc_iter)
#If MM THEORY (not just pointcharges)
if mm_theory is not None:
#Sanity check. Same number of atoms in fragment and MM object ?
if fragment.numatoms != mm_theory.numatoms:
print("")
print(BC.FAIL,"Number of atoms in fragment ({}) and MMtheory object differ ({})".format(fragment.numatoms,mm_theory.numatoms),BC.END)
print(BC.FAIL,"This does not make sense. Check coordinates and forcefield files. Exiting...", BC.END)
ashexit()
#Add possible exception for QM-QM atoms here.
#Maybe easier to just just set charges to 0. LJ for QM-QM still needs to be done by MM code
#if self.mm_theory_name == "OpenMMTheory":
#
#print("Now adding exceptions for frozen atoms")
#if len(self.frozenatoms) > 0:
# print("Here adding exceptions for OpenMM")
# print("Frozen-atom exceptions currently inactive...")
#print("Num frozen atoms: ", len(self.frozenatoms))
#Disabling for now, since so bloody slow. Need to speed up
#mm_theory.addexceptions(self.frozenatoms)
#Check if we need linkatoms by getting boundary atoms dict:
#Update: Tolerance modification to make sure we definitely catch connected atoms and get QM-MM boundary right.
#Scale=1.0 and tol=0.1 fails for S-C bond in rubredoxin from a classical MD run
#Bumping up a bit here.
#21 Sep 2023. bumping from +0.1 to +0.2. C-C bond in lysine failed
conn_scale=ash.settings_ash.settings_dict["scale"]
conn_tolerance=ash.settings_ash.settings_dict["tol"]+0.2
#If QM-MM boundary issue and ASH exits then printing QM-coordinates is useful
print("QM-region coordinates (before linkatoms):")
ash.modules.module_coords.print_coords_for_atoms(self.coords, self.elems, self.qmatoms, labels=self.qmatoms)
print()
self.boundaryatoms = ash.modules.module_coords.get_boundary_atoms(self.qmatoms, self.coords, self.elems, conn_scale,
conn_tolerance, excludeboundaryatomlist=excludeboundaryatomlist, unusualboundary=unusualboundary)
if len(self.boundaryatoms) >0:
print("Found covalent QM-MM boundary. Linkatoms option set to True")
print("Boundaryatoms (QM:MM pairs):", self.boundaryatoms)
print("Note: used connectivity settings, scale={} and tol={} to determine boundary.".format(conn_scale,conn_tolerance))
self.linkatoms=True
#Get MM boundary information. Stored as self.MMboundarydict
self.get_MMboundary(conn_scale,conn_tolerance)
else:
print("No covalent QM-MM boundary. Linkatoms option set to False")
self.linkatoms=False
#Removing possible QM atom constraints in OpenMMTheory
#Will only apply when running OpenMM_Opt or OpenMM_MD
if self.mm_theory_name == "OpenMMTheory":
self.mm_theory.remove_constraints_for_atoms(self.qmatoms)
#Remove bonded interactions in MM part. Only in OpenMM. Assuming they were never defined in NonbondedTHeory
print("Removing bonded terms for QM-region in MMtheory")
self.mm_theory.modify_bonded_forces(self.qmatoms)
#NOTE: Temporary. Adding exceptions for nonbonded QM atoms. Will ignore QM-QM Coulomb and LJ interactions.
#NOTE: For QM-MM interactions Coulomb charges are zeroed below (update_charges and delete_exceptions)
print("Removing nonbonded terms for QM-region in MMtheory (QM-QM interactions)")
self.mm_theory.addexceptions(self.qmatoms)
########################
#CHANGE CHARGES
########################
# Keeping self.charges as originally defined.
#Setting QM charges to 0 since electrostatic embedding
#and Charge-shift QM-MM boundary
#Zero QM charges for electrostatic embedding
#TODO: DO here or inside run instead?? Needed for MM code.
if self.embedding=="Elstat":
print("Charges of QM atoms set to 0 (since Electrostatic Embedding):")
self.ZeroQMCharges() #Modifies self.charges_qmregionzeroed
#print("length of self.charges_qmregionzeroed :", len(self.charges_qmregionzeroed))
#TODO: make sure this works for OpenMM and for NonBondedTheory
# Updating charges in MM object.
self.mm_theory.update_charges(self.qmatoms,[0.0 for i in self.qmatoms])
if self.mm_theory_name == "OpenMMTheory":
#Deleting Coulomb exception interactions involving QM and MM atoms
self.mm_theory.delete_exceptions(self.qmatoms)
#Option to create OpenMM externalforce that handles full system
if self.openmm_externalforce == True:
print("openmm_externalforce is True")
#print("Creating new OpenMM custom external force")
#MOVED FROM HERE TO OPENMM_MD
#Printing charges: all or only QM
if self.printlevel > 2:
for i in self.allatoms:
if i in self.qmatoms:
if self.embedding=="Elstat":
print("QM atom {} ({}) charge: {}".format(i, self.elems[i], self.charges_qmregionzeroed[i]))
else:
print("QM atom {} ({}) charge: {}".format(i, self.elems[i], self.charges[i]))
else:
if self.printlevel > 3:
print("MM atom {} ({}) charge: {}".format(i, self.elems[i], self.charges_qmregionzeroed[i]))
blankline()
else:
#Case: No actual MM theory but we still want to zero charges for QM elstate embedding calculation
#TODO: Remove option for no MM theory or keep this ??
if self.embedding=="Elstat":
self.ZeroQMCharges() #Modifies self.charges_qmregionzeroed
print_time_rel(module_init_time, modulename='QM/MM object creation', currprintlevel=self.printlevel, moduleindex=3)
#From QM1:MM1 boundary dict, get MM1:MMx boundary dict (atoms connected to MM1)
def get_MMboundary(self,scale,tol):
timeA=time.time()
# if boundarydict is not empty we need to zero MM1 charge and distribute charge from MM1 atom to MM2,MM3,MM4
#Creating dictionary for each MM1 atom and its connected atoms: MM2-4
self.MMboundarydict={}
for (QM1atom,MM1atom) in self.boundaryatoms.items():
connatoms = ash.modules.module_coords.get_connected_atoms(self.coords, self.elems, scale,tol, MM1atom)
#Deleting QM-atom from connatoms list
connatoms.remove(QM1atom)
self.MMboundarydict[MM1atom] = connatoms
print("")
print("MM boundary (MM1:MMx pairs):", self.MMboundarydict)
print_time_rel(timeA, modulename="get_MMboundary")
# Set QMcharges to Zero and shift charges at boundary
#TODO: Add both L2 scheme (delete whole charge-group of M1) and charge-shifting scheme (shift charges to Mx atoms and add dipoles for each Mx atom)
def ZeroQMCharges(self):
timeA=time.time()
print("Setting QM charges to Zero")
#Looping over charges and setting QM atoms to zero
#1. Copy charges to charges_qmregionzeroed
self.charges_qmregionzeroed=copy.copy(self.charges)
#2. change charge for QM-atom
for i, c in enumerate(self.charges_qmregionzeroed):
#Setting QMatom charge to 0
if i in self.qmatoms:
self.charges_qmregionzeroed[i] = 0.0
#3. Flag that this has been done
self.QMChargesZeroed = True
print_time_rel(timeA, modulename="ZeroQMCharges")
def ShiftMMCharges(self):
timeA=time.time()
if self.printlevel > 1:
print("Shifting MM charges at QM-MM boundary.")
#print("len self.charges_qmregionzeroed: ", len(self.charges_qmregionzeroed))
#print("len self.charges: ", len(self.charges))
#Create self.pointcharges list
self.pointcharges=copy.copy(self.charges_qmregionzeroed)
#Looping over charges and setting QM/MM1 atoms to zero and shifting charge to neighbouring atoms
for i, c in enumerate(self.pointcharges):
#If index corresponds to MMatom at boundary, set charge to 0 (charge-shifting
if i in self.MMboundarydict.keys():
MM1charge = self.charges[i]
#print("MM1atom charge: ", MM1charge)
self.pointcharges[i] = 0.0
#MM1 charge fraction to be divided onto the other MM atoms
MM1charge_fract = MM1charge / len(self.MMboundarydict[i])
#print("MM1charge_fract :", MM1charge_fract)
#Putting the fractional charge on each MM2 atom
for MMx in self.MMboundarydict[i]:
#print("MMx : ", MMx)
#print("Old charge : ", self.charges_qmregionzeroed[MMx])
self.pointcharges[MMx] += MM1charge_fract
#print("New charge : ", self.charges_qmregionzeroed[MMx])
#exit()
#print_time_rel(timeA, modulename="ShiftMMCharges", currprintlevel=self.printlevel, currthreshold=1)
#Create dipole charge (twice) for each MM2 atom that gets fraction of MM1 charge
def get_dipole_charge(self,delq,direction,mm1index,mm2index):
#timeA=time.time()
#Distance between MM1 and MM2
MM_distance = ash.modules.module_coords.distance_between_atoms(fragment=self.fragment, atoms=[mm1index, mm2index])
#Coordinates
mm1coords=np.array(self.fragment.coords[mm1index])
mm2coords=np.array(self.fragment.coords[mm2index])
SHIFT=0.15
#Normalize vector
def vnorm(p1):
r = math.sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]))
v1=np.array([p1[0] / r, p1[1] / r, p1[2] /r])
return v1
diffvector=mm2coords-mm1coords
normdiffvector=vnorm(diffvector)
#Dipole
d = delq*2.5
#Charge (abs value)
q0 = 0.5 * d / SHIFT
#print("q0 : ", q0)
#Actual shift
#print("direction : ", direction)
shift = direction * SHIFT * ( MM_distance / 2.5 )
#print("shift : ", shift)
#Position
#print("normdiffvector :", normdiffvector)
#print(normdiffvector*shift)
pos = mm2coords+np.array((shift*normdiffvector))
#print("pos :", pos)
#Returning charge with sign based on direction and position
#Return coords as regular list
#print_time_rel(timeA, modulename="get_dipole_charge")
return -q0*direction,list(pos)
def SetDipoleCharges(self):
timeA=time.time()
if self.printlevel > 1:
print("Adding extra charges to preserve dipole moment for charge-shifting")
#Adding 2 dipole pointcharges for each MM2 atom
self.dipole_charges = []
self.dipole_coords = []
for MM1,MMx in self.MMboundarydict.items():
#Getting original MM1 charge (before set to 0)
MM1charge = self.charges[MM1]
MM1charge_fract=MM1charge/len(MMx)
for MM in MMx:
q_d1, pos_d1 = self.get_dipole_charge(MM1charge_fract,1,MM1,MM)
q_d2, pos_d2 = self.get_dipole_charge(MM1charge_fract,-1,MM1,MM)
self.dipole_charges.append(q_d1)
self.dipole_charges.append(q_d2)
self.dipole_coords.append(pos_d1)
self.dipole_coords.append(pos_d2)
#print_time_rel(timeA, modulename="SetDipoleCharges", currprintlevel=self.printlevel, currthreshold=1)
#More efficient version than previous loop
def make_QM_PC_gradient(self):
self.QM_PC_gradient = np.zeros((len(self.allatoms), 3))
qmatom_indices = np.where(np.isin(self.allatoms, self.qmatoms))[0]
pcatom_indices = np.where(~np.isin(self.allatoms, self.qmatoms))[0] # ~ is NOT operator in numpy
self.QM_PC_gradient[qmatom_indices] = self.QMgradient_wo_linkatoms[:len(qmatom_indices)]
self.QM_PC_gradient[pcatom_indices] = self.PCgradient[:len(pcatom_indices)]
return
#TruncatedPCfunction control flow for pointcharge field passed to QM program
def TruncatedPCfunction(self):
self.TruncatedPCcalls+=1
print("TruncatedPC approximation!")
if self.TruncatedPCcalls == 1:
self.TruncatedPC_recalc_flag=True
print("This is first QM/MM run. Will calculate Full-Trunc correction in this step")
#Origin coords point is center of QM-region
origincoords=ash.modules.module_coords.get_centroid(self.qmcoords)
#Determine the indices associated with the truncated PC field once
self.determine_truncatedPC_indices(origincoords)
print("Truncated PC-region size: {} charges".format(len(self.truncated_PC_region_indices)))
#Saving full PCs and coords for 1st iteration
self.pointcharges_full=copy.copy(self.pointcharges)
self.pointchargecoords_full=copy.copy(self.pointchargecoords)
#Determining truncated PC-field
self.pointcharges=[self.pointcharges[i] for i in self.truncated_PC_region_indices]
self.pointchargecoords=np.take(self.pointchargecoords, self.truncated_PC_region_indices, axis=0)
elif self.TruncatedPCcalls % self.TruncatedPC_recalc_iter == 0:
self.TruncatedPC_recalc_flag=True
print(f"This is QM/MM run no. {self.TruncatedPCcalls}. Will calculate Full-Trunc correction in this step")
print("self.TruncatedPC_recalc_iter:", self.TruncatedPC_recalc_iter)
origincoords=ash.modules.module_coords.get_centroid(self.qmcoords)
#Determine the indices associated with the truncated PC field once
self.determine_truncatedPC_indices(origincoords)
print("Truncated PC-region size: {} charges".format(len(self.truncated_PC_region_indices)))
#Saving full PCs and coords for 1st iteration
self.pointcharges_full=copy.copy(self.pointcharges)
self.pointchargecoords_full=copy.copy(self.pointchargecoords)
#Determining truncated PC-field
self.pointcharges=[self.pointcharges[i] for i in self.truncated_PC_region_indices]
self.pointchargecoords=np.take(self.pointchargecoords, self.truncated_PC_region_indices, axis=0)
else:
self.TruncatedPC_recalc_flag=False
print("This is QM/MM run no. {}. Using approximate truncated PC field: {} charges".format(self.TruncatedPCcalls,len(self.truncated_PC_region_indices)))
self.pointcharges=[self.pointcharges[i] for i in self.truncated_PC_region_indices]
self.pointchargecoords=np.take(self.pointchargecoords, self.truncated_PC_region_indices, axis=0)
#Determine truncated PC field indices based on initial coordinates
#Coordinates and charges for each Opt cycle defined later.
def determine_truncatedPC_indices(self,origincoords):
region_indices=[]
for index,allc in enumerate(self.pointchargecoords):
dist=ash.modules.module_coords.distance(origincoords,allc)
if dist < self.TruncPCRadius:
region_indices.append(index)
#Only unique and sorting:
self.truncated_PC_region_indices = np.unique(region_indices).tolist()
#Removing dipole charges also (end of list)
#NOTE:
#num_dipole_charges=len(self.dipole_charges)
#self.truncated_PC_region_indices = self.truncated_PC_region_indices[0:-num_dipole_charges]
def oldcalculate_truncPC_gradient_correction(self,QMgradient_full, PCgradient_full, QMgradient_trunc, PCgradient_trunc):
#Correction for QM-atom gradient
self.original_QMcorrection_gradient=np.zeros((len(QMgradient_full)-self.num_linkatoms, 3))
#Correction for PC gradient
self.original_PCcorrection_gradient=np.zeros((len(PCgradient_full), 3))
for i in range(0,len(QMgradient_full)-self.num_linkatoms):
#print("QM index:", i)
qmfullgrad=QMgradient_full[i]
qmtruncgrad=QMgradient_trunc[i]
#print("qmfullgrad:", qmfullgrad)
#print("qmtruncgrad:", qmtruncgrad)
qmdifference=qmfullgrad-qmtruncgrad
#print("qmdifference:", qmdifference)
self.original_QMcorrection_gradient[i] = qmdifference
count=0
for i in range(0,len(PCgradient_full)):
if i in self.truncated_PC_region_indices:
#print("TruncPC index:", i)
pcfullgrad=PCgradient_full[i]
pctruncgrad=PCgradient_trunc[count]
#print("pcfullgrad:", pcfullgrad)
#print("pctruncgrad:", pctruncgrad)
difference=pcfullgrad-pctruncgrad
#print("difference:", difference)
self.original_PCcorrection_gradient[i] = difference
count+=1
else:
# Keep original full contribution
self.original_PCcorrection_gradient[i] = PCgradient_full[i]
return
#New more efficient version
def calculate_truncPC_gradient_correction(self, QMgradient_full, PCgradient_full, QMgradient_trunc, PCgradient_trunc):
#QM part
qm_difference = QMgradient_full[:len(QMgradient_full)-self.num_linkatoms] - QMgradient_trunc[:len(QMgradient_full)-self.num_linkatoms]
self.original_QMcorrection_gradient = qm_difference
#PC part
truncated_indices = np.array(self.truncated_PC_region_indices)
pc_difference = np.zeros((len(PCgradient_full), 3))
pc_difference[truncated_indices] = PCgradient_full[truncated_indices] - PCgradient_trunc
pc_difference[~np.isin(np.arange(len(PCgradient_full)), truncated_indices)] = PCgradient_full[~np.isin(np.arange(len(PCgradient_full)), truncated_indices)]
self.original_PCcorrection_gradient = pc_difference
return
#This updates the calculated truncated PC gradient to be full-system gradient
#by combining with the original 1st step correction
def oldTruncatedPCgradientupdate(self,QMgradient_wo_linkatoms,PCgradient):
#QM part
newQMgradient_wo_linkatoms = QMgradient_wo_linkatoms + self.original_QMcorrection_gradient
#PC part
new_full_PC_gradient=np.zeros((len(self.original_PCcorrection_gradient), 3))
count=0
for i in range(0,len(new_full_PC_gradient)):
if i in self.truncated_PC_region_indices:
#Now updating with gradient from active region
new_full_PC_gradient[i] = self.original_PCcorrection_gradient[i] + PCgradient[count]
count+=1
else:
new_full_PC_gradient[i] = self.original_PCcorrection_gradient[i]
return newQMgradient_wo_linkatoms, new_full_PC_gradient
def TruncatedPCgradientupdate(self, QMgradient_wo_linkatoms, PCgradient):
newQMgradient_wo_linkatoms = QMgradient_wo_linkatoms + self.original_QMcorrection_gradient
new_full_PC_gradient = np.copy(self.original_PCcorrection_gradient)
new_full_PC_gradient[self.truncated_PC_region_indices] += PCgradient
return newQMgradient_wo_linkatoms, new_full_PC_gradient
def set_numcores(self,numcores):
print(f"Setting new numcores {numcores}for QMtheory and MMtheory")
self.qm_theory.set_numcores(numcores)
self.mm_theory.set_numcores(numcores)
#Method to grab dipole moment from outputfile (assumes run has been executed)
def get_dipole_moment(self):
try:
print("Grabbing dipole moment from QM-part of QM/MM theory.")
dipole = self.qm_theory.get_dipole_moment()
except:
print("Error: Could not grab dipole moment from QM-part of QM/MM theory.")
return dipole
#Method to polarizability from outputfile (assumes run has been executed)
def get_polarizability_tensor(self):
try:
print("Grabbing polarizability from QM-part of QM/MM theory.")
polarizability = self.qm_theory.get_polarizability_tensor()
except:
print("Error: Could not grab polarizability from QM-part of QM/MM theory.")
return polarizability
#General run
def run(self, current_coords=None, elems=None, Grad=False, numcores=1, exit_after_customexternalforce_update=False, label=None, charge=None, mult=None):
if self.embedding == "Mechanical":
return self.mech_run(current_coords=current_coords, elems=elems, Grad=Grad, numcores=numcores, exit_after_customexternalforce_update=exit_after_customexternalforce_update,
label=label, charge=charge, mult=mult)
elif self.embedding == "Elstat":
return self.elstat_run(current_coords=current_coords, elems=elems, Grad=Grad, numcores=numcores, exit_after_customexternalforce_update=exit_after_customexternalforce_update,
label=label, charge=charge, mult=mult)
else:
print("Unknown embedding. Exiting")
ashexit()
#Mechanical embedding run
def mech_run(self, current_coords=None, elems=None, Grad=False, numcores=1, exit_after_customexternalforce_update=False, label=None, charge=None, mult=None):
print("Mechanical embedding not ready yet")
#Since elstat-run is so complicated it is easier to just write new mechanical run
#mechanical embedding may also come with its own complexities
#NOTE: We should also reduce complexity of elstat-run by moving code into methods
ashexit()
#Electrostatic embedding run
def elstat_run(self, current_coords=None, elems=None, Grad=False, numcores=1, exit_after_customexternalforce_update=False, label=None, charge=None, mult=None):
module_init_time=time.time()
CheckpointTime = time.time()
if self.printlevel >= 2:
print(BC.WARNING, BC.BOLD, "------------RUNNING QM/MM MODULE-------------", BC.END)
print("QM Module:", self.qm_theory_name)
print("MM Module:", self.mm_theory_name)
#OPTION: QM-region charge/mult from QMMMTheory definition
#If qm_charge/qm_mult defined then we use. Otherwise charge/mult may have been defined by jobtype-function and passed on via run
if self.qm_charge != None:
if self.printlevel > 1: print("Charge provided from QMMMTheory object: ", self.qm_charge)
charge=self.qm_charge
if self.qm_mult != None:
if self.printlevel > 1: print("Mult provided from QMMMTheory object: ", self.qm_mult)
mult=self.qm_mult
#Checking if charge and mult has been provided. Exit if not.
if charge == None or mult == None:
print(BC.FAIL, "Error. charge and mult has not been defined for QMMMTheory.run method", BC.END)
ashexit()
if self.printlevel >1 : print("QM-region Charge: {} Mult: {}".format(charge,mult))
#If no coords provided to run (from Optimizer or NumFreq or MD) then use coords associated with object.
#if len(current_coords) != 0:
if current_coords is not None:
pass
else:
current_coords=self.coords
if self.embedding=="Elstat":
PC=True
else:
PC=False
#If numcores was set when calling .run then using, otherwise use self.numcores
if numcores==1:
numcores=self.numcores
if self.printlevel >= 2:
print("Running QM/MM object with {} cores available".format(numcores))
#Updating QM coords and MM coords.
#TODO: Should we use different name for updated QMcoords and MMcoords here??
#self.qmcoords=[current_coords[i] for i in self.qmatoms]
#self.mmcoords=[current_coords[i] for i in self.mmatoms]
self.qmcoords=np.take(current_coords,self.qmatoms,axis=0)
self.mmcoords=np.take(current_coords,self.mmatoms,axis=0)
self.qmelems=[self.elems[i] for i in self.qmatoms]
self.mmelems=[self.elems[i] for i in self.mmatoms]
#LINKATOMS
#1. Get linkatoms coordinates
if self.linkatoms==True:
linkatoms_dict = ash.modules.module_coords.get_linkatom_positions(self.boundaryatoms,self.qmatoms, current_coords, self.elems)
printdebug("linkatoms_dict:", linkatoms_dict)
#2. Add linkatom coordinates to qmcoords???
if self.printlevel > 1: print("Adding linkatom positions to QM coords")
linkatoms_indices=[]
#Sort by QM atoms:
printdebug("linkatoms_dict.keys :", linkatoms_dict.keys())
for pair in sorted(linkatoms_dict.keys()):
printdebug("Pair :", pair)
#self.qmcoords.append(linkatoms_dict[pair])
self.qmcoords = np.append(self.qmcoords,np.array([linkatoms_dict[pair]]), axis=0)
#print("self.qmcoords :", self.qmcoords)
#print(len(self.qmcoords))
#ashexit()
#Linkatom indices for book-keeping
linkatoms_indices.append(len(self.qmcoords)-1)
printdebug("linkatoms_indices: ", linkatoms_indices)
self.num_linkatoms=len(linkatoms_indices)
if self.printlevel > 1: print(f"There are {self.num_linkatoms} linkatoms")
#TODO: Modify qm_elems list. Use self.qmelems or separate qmelems ?
#TODO: Should we do this at object creation instead?
current_qmelems=self.qmelems + ['H']*len(linkatoms_dict)
#print("")
#print("current_qmelems :", current_qmelems)
#Do Charge-shifting. MM1 charge distributed to MM2 atoms
if self.embedding=="Elstat":
if self.printlevel > 1: print("Doing charge-shifting...")
#print("Before: self.pointcharges are: ", self.pointcharges)
self.ShiftMMCharges() # Creates self.pointcharges
#print("After: self.pointcharges are: ", self.pointcharges)
if self.printlevel > 1: print("Number of pointcharges defined for whole system: ", len(self.pointcharges))
#TODO: Code alternative to Charge-shifting: L2 scheme which deletes whole charge-group that MM1 belongs to
# Defining pointcharges as only containing MM atoms
if self.printlevel > 1: print("Number of MM atoms:", len(self.mmatoms))
self.pointcharges=[self.pointcharges[i] for i in self.mmatoms]
#print("After: self.pointcharges are: ", self.pointcharges)
if self.printlevel > 1: print("Number of pointcharges defined for MM region: ", len(self.pointcharges))
#Set
if self.dipole_correction is True:
print("Dipole correction is on. Adding dipole charges")
self.SetDipoleCharges() #Creates self.dipole_charges and self.dipole_coords
#Adding dipole charge coords to MM coords (given to QM code) and defining pointchargecoords
if self.printlevel > 1: print("Adding {} dipole charges to PC environment".format(len(self.dipole_charges)))
self.pointchargecoords=np.append(self.mmcoords,np.array(self.dipole_coords), axis=0)
#Adding dipole charges to MM charges list (given to QM code)
self.pointcharges=self.pointcharges+self.dipole_charges
#Using element H for dipole charges. Only matters for CP2K GEEP
mm_elems_for_qmprogram=self.mmelems + ['H']*len(self.dipole_charges)
if self.printlevel > 1: print("Number of pointcharges after dipole addition: ", len(self.pointcharges))
else:
print("Dipole correction is off. Not adding any dipole charges")
self.pointchargecoords = self.mmcoords
mm_elems_for_qmprogram=self.mmelems
if self.printlevel > 1: print("Number of pointcharges: ", len(self.pointcharges))
else:
mm_elems_for_qmprogram=self.mmelems
self.num_linkatoms=0
#If no linkatoms then use original self.qmelems
current_qmelems = self.qmelems
#If no linkatoms then self.pointcharges are just original charges with QM-region zeroed
if self.embedding=="Elstat":
self.pointcharges=[self.charges_qmregionzeroed[i] for i in self.mmatoms]
#If no linkatoms MM coordinates are the same
self.pointchargecoords=self.mmcoords
#NOTE: Now we have updated MM-coordinates (if doing linkatoms, wtih dipolecharges etc) and updated mm-charges (more, due to dipolecharges if linkatoms)
# We also have MMcharges that have been set to zero due to QM/mm
#We do not delete charges but set to zero
if self.printlevel > 1:
print("Number of charges:", len(self.pointcharges))
print("Number of charge coordinates:", len(self.pointchargecoords))
#TRUNCATED PC Option: Speeding up QM/MM jobs of large systems by passing only a truncated PC field to the QM-code most of the time
# Speeds up QM-pointcharge gradient that otherwise dominates
if self.TruncatedPC is True:
self.TruncatedPCfunction()
#Modifies self.pointcharges and self.pointchargecoords
#print("Number of charges after truncation :", len(self.pointcharges))
#print("Number of charge coordinates after truncation :", len(self.pointchargecoords))
#If no qmatoms then do MM-only
if len(self.qmatoms) == 0:
print("No qmatoms list provided. Setting QMtheory to None")
self.qm_theory_name="None"
self.QMenergy=0.0
#print_time_rel(CheckpointTime, modulename='QM/MM run prep', moduleindex=2, currprintlevel=self.printlevel, currthreshold=1)
if self.printlevel > 1: print("Number of PCs provided to QM-program:", len(self.pointcharges))
#QM-part
if self.qm_theory_name == "None" or self.qm_theory_name == "ZeroTheory":
print("No QMtheory. Skipping QM calc")
QMenergy=0.0;self.linkatoms=False;PCgradient=np.array([0.0, 0.0, 0.0])
QMgradient=np.array([0.0, 0.0, 0.0])
#General QM-code energy+gradient call.
#TODO: Add check whether QM-code supports both pointcharges and pointcharge-gradient?
else:
#Calling QM theory, providing current QM and MM coordinates.
if Grad==True:
if PC==True:
QMenergy, QMgradient, PCgradient = self.qm_theory.run(current_coords=self.qmcoords,
current_MM_coords=self.pointchargecoords,
MMcharges=self.pointcharges,
qm_elems=current_qmelems, mm_elems=mm_elems_for_qmprogram,
charge=charge, mult=mult,
Grad=True, PC=True, numcores=numcores)
#shutil.copyfile(self.qm_theory.filename+'.out', self.qm_theory.filename+'_full'+'.out')
else:
QMenergy, QMgradient = self.qm_theory.run(current_coords=self.qmcoords,
current_MM_coords=self.pointchargecoords, MMcharges=self.pointcharges,
qm_elems=current_qmelems, Grad=True, PC=False, numcores=numcores, charge=charge, mult=mult)
else:
QMenergy = self.qm_theory.run(current_coords=self.qmcoords,
current_MM_coords=self.pointchargecoords, MMcharges=self.pointcharges,
qm_elems=current_qmelems, Grad=False, PC=PC, numcores=numcores, charge=charge, mult=mult)
print_time_rel(CheckpointTime, modulename='QM step', moduleindex=2,currprintlevel=self.printlevel, currthreshold=1)
CheckpointTime = time.time()
#Final QM/MM gradient. Combine QM gradient, MM gradient, PC-gradient (elstat MM gradient from QM code).
# Do linkatom force projections in the end
#UPDATE: Do MM step in the end now so that we have options for OpenMM extern force
if Grad == True:
Grad_prep_CheckpointTime = time.time()
#assert len(self.allatoms) == len(self.MMgradient)
#Defining QMgradient without linkatoms if present
if self.linkatoms==True:
self.QMgradient = QMgradient
QMgradient_wo_linkatoms=QMgradient[0:-self.num_linkatoms] #remove linkatoms
else:
self.QMgradient = QMgradient
QMgradient_wo_linkatoms=QMgradient
#if self.printlevel >= 2:
# ash.modules.module_coords.write_coords_all(self.QMgradient_wo_linkatoms, self.qmelems, indices=self.allatoms, file="QMgradient_wo_linkatoms", description="QM+ gradient withoutlinkatoms (au/Bohr):")
#TRUNCATED PC Option:
if self.TruncatedPC is True:
#DONE ONCE: CALCULATE FULL PC GRADIENT TO DETERMINE CORRECTION
if self.TruncatedPC_recalc_flag is True:
CheckpointTime = time.time()
truncfullCheckpointTime = time.time()
#We have calculated truncated QM and PC gradient
QMgradient_trunc = QMgradient
PCgradient_trunc = PCgradient
print("Now calculating full QM and PC gradient")
print("Number of PCs provided to QM-program:", len(self.pointcharges_full))
QMenergy_full, QMgradient_full, PCgradient_full = self.qm_theory.run(current_coords=self.qmcoords,
current_MM_coords=self.pointchargecoords_full,
MMcharges=self.pointcharges_full,
qm_elems=current_qmelems, charge=charge, mult=mult,
Grad=True, PC=True, numcores=numcores)
print_time_rel(CheckpointTime, modulename='trunc-pc full calculation', moduleindex=3)
CheckpointTime = time.time()
#try:
# shutil.copyfile(self.qm_theory.filename+'.out', self.qm_theory.filename+'_full'+'.out')
#except:
# pass
#TruncPC correction to QM energy
self.truncPC_E_correction = QMenergy_full - QMenergy
print(f"Truncated PC energy correction: {self.truncPC_E_correction} Eh")
self.QMenergy = QMenergy + self.truncPC_E_correction
#Now determine the correction once and for all
CheckpointTime = time.time()
self.calculate_truncPC_gradient_correction(QMgradient_full, PCgradient_full, QMgradient_trunc, PCgradient_trunc)
print_time_rel(CheckpointTime, modulename='calculate_truncPC_gradient_correction', moduleindex=3)
CheckpointTime = time.time()
#Now defining final QMgradient and PCgradient
self.QMgradient_wo_linkatoms, self.PCgradient = self.TruncatedPCgradientupdate(QMgradient_wo_linkatoms,PCgradient)
print_time_rel(CheckpointTime, modulename='truncPC_gradient update ', moduleindex=3)
print_time_rel(truncfullCheckpointTime, modulename='trunc-full-step pcgrad update', moduleindex=3)
else:
CheckpointTime = time.time()
#TruncPC correction to QM energy
self.QMenergy = QMenergy + self.truncPC_E_correction
self.QMgradient_wo_linkatoms, self.PCgradient = self.TruncatedPCgradientupdate(QMgradient_wo_linkatoms,PCgradient)
print_time_rel(CheckpointTime, modulename='trunc pcgrad update', moduleindex=3)
else:
self.QMenergy = QMenergy
#No TruncPC approximation active. No change to original QM and PCgradient from QMcode
self.QMgradient_wo_linkatoms = QMgradient_wo_linkatoms
if self.embedding=="Elstat":
self.PCgradient = PCgradient
#Initialize QM_PC gradient (has full system size) and fill up
#TODO: This can be made more efficient
CheckpointTime = time.time()
self.make_QM_PC_gradient() #populates self.QM_PC_gradient
print_time_rel(CheckpointTime, modulename='QMpcgrad prepare', moduleindex=3, currprintlevel=self.printlevel, currthreshold=2)
#self.QM_PC_gradient = np.zeros((len(self.allatoms), 3))
#qmcount=0;pccount=0
#for i in self.allatoms:
# if i in self.qmatoms:
# #QM-gradient. Linkatom gradients are skipped
# self.QM_PC_gradient[i]=self.QMgradient_wo_linkatoms[qmcount]
# qmcount+=1
# else:
# #Pointcharge-gradient. Dipole-charge gradients are skipped (never reached)
# self.QM_PC_gradient[i] = self.PCgradient[pccount]
# pccount += 1
#assert qmcount == len(self.qmatoms)
#assert pccount == len(self.mmatoms)
#assert qmcount+pccount == len(self.allatoms)
#print(" qmcount+pccount:", qmcount+pccount)
#print("len(self.allatoms):", len(self.allatoms))
#print("len self.QM_PC_gradient", len(self.QM_PC_gradient))
#ash.modules.module_coords.write_coords_all(self.QM_PC_gradient, self.elems, indices=self.allatoms, file="QM_PC_gradient", description="QM_PC_gradient (au/Bohr):")
#if self.printlevel >= 2:
# ash.modules.module_coords.write_coords_all(self.PCgradient, self.mmatoms, indices=self.allatoms, file="PCgradient", description="PC gradient (au/Bohr):")
#if self.printlevel >= 2:
# ash.modules.module_coords.write_coords_all(self.QM_PC_gradient, self.elems, indices=self.allatoms, file="QM+PCgradient_before_linkatomproj", description="QM+PC gradient before linkatomproj (au/Bohr):")
#LINKATOM FORCE PROJECTION
# Add contribution to QM1 and MM1 contribution???
if self.linkatoms==True:
CheckpointTime = time.time()
#print("here")
#print("linkatoms_dict: ", linkatoms_dict)
#print("linkatoms_indices: ", linkatoms_indices)
for pair in sorted(linkatoms_dict.keys()):
printdebug("pair: ", pair)
#Grabbing linkatom data
linkatomindex=linkatoms_indices.pop(0)
printdebug("linkatomindex:", linkatomindex)
Lgrad=self.QMgradient[linkatomindex]
printdebug("Lgrad:",Lgrad)
Lcoord=linkatoms_dict[pair]
printdebug("Lcoord:", Lcoord)
#Grabbing QMatom info
fullatomindex_qm=pair[0]
printdebug("fullatomindex_qm:", fullatomindex_qm)
printdebug("self.qmatoms:", self.qmatoms)
qmatomindex=fullindex_to_qmindex(fullatomindex_qm,self.qmatoms)
printdebug("qmatomindex:", qmatomindex)
Qcoord=self.qmcoords[qmatomindex]
printdebug("Qcoords: ", Qcoord)
Qgrad=self.QM_PC_gradient[fullatomindex_qm]
printdebug("Qgrad (full QM/MM grad)s:", Qgrad)
#Grabbing MMatom info
fullatomindex_mm=pair[1]
printdebug("fullatomindex_mm:", fullatomindex_mm)
Mcoord=current_coords[fullatomindex_mm]
printdebug("Mcoord:", Mcoord)
Mgrad=self.QM_PC_gradient[fullatomindex_mm]
printdebug("Mgrad (full QM/MM grad): ", Mgrad)
#Now grabbed all components, calculating new projected gradient on QM atom and MM atom
newQgrad,newMgrad = linkatom_force_fix(Qcoord, Mcoord, Lcoord, Qgrad,Mgrad,Lgrad)
printdebug("newQgrad: ", newQgrad)
printdebug("newMgrad: ", newMgrad)
#Updating full QM_PC_gradient (used to be QM/MM gradient)
#self.QM_MM_gradient[fullatomindex_qm] = newQgrad
#self.QM_MM_gradient[fullatomindex_mm] = newMgrad
self.QM_PC_gradient[fullatomindex_qm] = newQgrad
self.QM_PC_gradient[fullatomindex_mm] = newMgrad
print_time_rel(CheckpointTime, modulename='linkatomgrad prepare', moduleindex=2, currprintlevel=self.printlevel, currthreshold=1)
print_time_rel(Grad_prep_CheckpointTime, modulename='QM/MM gradient prepare', moduleindex=2, currprintlevel=self.printlevel, currthreshold=1)
CheckpointTime = time.time()
else:
#No Grad
self.QMenergy = QMenergy
#if self.printlevel >= 2:
# ash.modules.module_coords.write_coords_all(self.QM_PC_gradient, self.elems, indices=self.allatoms, file="QM+PCgradient", description="QM+PC gradient (au/Bohr):")
# MM THEORY
if self.mm_theory_name == "NonBondedTheory":
if self.printlevel >= 2:
print("Running MM theory as part of QM/MM.")
print("Using MM on full system. Charges for QM region have to be set to zero ")
#printdebug("Charges for full system is: ", self.charges)
print("Passing QM atoms to MMtheory run so that QM-QM pairs are skipped in pairlist")
print("Passing active atoms to MMtheory run so that frozen pairs are skipped in pairlist")
assert len(current_coords) == len(self.charges_qmregionzeroed)
# NOTE: charges_qmregionzeroed for full system but with QM-charges zeroed (no other modifications)
#NOTE: Using original system coords here (not with linkatoms, dipole etc.). Also not with deleted zero-charge coordinates.
#charges list for full system, can be zeroed but we still want the LJ interaction
self.MMenergy, self.MMgradient= self.mm_theory.run(current_coords=current_coords,
charges=self.charges_qmregionzeroed, connectivity=self.connectivity,
qmatoms=self.qmatoms, actatoms=self.actatoms)
elif self.mm_theory_name == "OpenMMTheory":
if self.printlevel >= 2:
print("Using OpenMM theory as part of QM/MM.")
if self.QMChargesZeroed==True:
if self.printlevel >= 2:
print("Using MM on full system. Charges for QM region {} have been set to zero ".format(self.qmatoms))
else:
print("QMCharges have not been zeroed")
ashexit()
#printdebug("Charges for full system is: ", self.charges)