Dear @chengsu ,
sorry this is a bug in the code. The code silently discards the “Hmag” keyword argument.
I’ve created an issue for this:
opened 01:14PM - 19 Aug 24 UTC
Bug
**Original reporter:** chengsu [via Forum](https://forum.mantidproject.org/t/que… stions-on-magnetic-moment-calculation-using-crystalfield/898)
**Describe the bug**
The Crystal Field physical properties calculations using the `getMagneticMoment` method for the `M(T)` case (`typeid==4`) is incorrect because it [discards the `Hmag` field](https://github.com/mantidproject/mantid/blob/78d0301522e6a4019ac7213e7a82c2fda96f1126/scripts/Inelastic/CrystalField/function.py#L705) - this line should be something like:
```python
if self._typeid == self.MAGNETISATION or self._typeid == self.MAGNETICMOMENT:
out += (",Hmag=%s" % (self._hmag))
```
This bug means that whatever value the user puts in for the `Hmag` keyword argument, they will obtain the same calculated curve, which is incorrect.
**To Reproduce**
```python
from mantid.simpleapi import CreateWorkspace
from CrystalField import CrystalField, PhysicalProperties
import numpy as np
cf = CrystalField('Ce', 'C4', Temperature=1, B20=0.2, B40=0.01, B44=-0.005)
T = np.linspace(1,300,100)
moment_T0p1 = CreateWorkspace(*cf.getMagneticMoment (0.1, Temperature = T , Hdir=[1,0,0], Hmag=0.1, Unit='bohr'))
moment_T1 = CreateWorkspace(*cf.getMagneticMoment (1, Temperature = T , Hdir=[1,0,0], Hmag=1, Unit='bohr'))
moment_T10 = CreateWorkspace(*cf.getMagneticMoment (100, Temperature = T , Hdir=[1,0,0], Hmag=10, Unit='bohr'))
```
Plot the three curves computed - they all lie on top of each other and they should not.
<!--Steps to reproduce the behavior:
For example
1. Go to '...'
2. Click on '....'
3. Scroll down to '....'
4. See error
-->
**Expected behavior**
The expected behaviour can be obtained using this work-around
```python
from mantid.simpleapi import CreateWorkspace
from CrystalField import CrystalField, PhysicalProperties
import numpy as np
def mymagmom(cf, **kwargs):
Hmag = kwargs.pop('Hmag')
if not isinstance(T, str) or 'mantid' in type(T):
twksp = CreateWorkspace(kwargs['Temperature'], kwargs['Temperature'], OutputWorkspace='_twksp')
else:
twksp = kwargs['Temperature']
funstr = cf.makePhysicalPropertiesFunction(PhysicalProperties("M(T)", **kwargs)) + f',Hmag={Hmag}'
return cf._calcSpectrum(funstr, twksp, 0)
moment_T0p1 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=0.1, Unit='bohr'))
moment_T1 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=1, Unit='bohr'))
moment_T10 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=10, Unit='bohr'))
```
You can now see that the three spectra do not overlap and that the original calculation corresponds to `Hmag=1`.
**Screenshots**
**Platform/Version (please complete the following information):**
- OS: All
- OS Version:
- Mantid Version: At least since 6.9
**Additional context**
so hopefully it will be fixed by the next release.
However, if you want to have something working now, you can use this work-around function:
from mantid.simpleapi import CreateWorkspace
from CrystalField import CrystalField, PhysicalProperties
import numpy as np
def mymagmom(cf, **kwargs):
Hmag = kwargs.pop('Hmag')
if not isinstance(T, str) or 'mantid' in type(T):
twksp = CreateWorkspace(kwargs['Temperature'], kwargs['Temperature'], OutputWorkspace='_twksp')
else:
twksp = kwargs['Temperature']
funstr = cf.makePhysicalPropertiesFunction(PhysicalProperties("M(T)", **kwargs)) + f',Hmag={Hmag}'
return cf._calcSpectrum(funstr, twksp, 0)
cf = CrystalField('Ce', 'C4', Temperature=1, B20=0.2, B40=0.01, B44=-0.005)
T = np.linspace(1,300,100)
moment2_T0p1 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=0.1, Unit='bohr'))
moment2_T1 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=1, Unit='bohr'))
moment2_T10 = CreateWorkspace(*mymagmom(cf, Temperature = T , Hdir=[1,0,0], Hmag=10, Unit='bohr'))