mulfit

Mulfit

The current version is 2.0.
The files are in /home/bio/Mulfit/mulfit-2.0
Documentations are in /home/bio/Mulfit/mulfit-2.0/data/mulfit.doc.

The program is compiled by running src/compile.sh - it worked well
The executable src/mulfit-2.0.exe is linked to /home/bio/Bin/Linux/mulfit

There is a script mulfit-2.0.sh showing how to run the program:
ln -s $1.mlt for014
ln -s $1.par for016
mulfit-2.0.exe
rm for014
rm for016
mv for017 $1.eml
mv for018 $1.sum

Thus, there are four files involved:
for014 (input) contains the multipoles
for016 (input) contains the input parameters (sample fie)
for017 (output) is the log file
for018 (output) is a summary of the results

The program can read two file formats: MULPOL = AMPAC, if the file (for014) starts with the line
DMASEM
and CADPAC otherwise. The former format is read by subroutine ramp, the latter by rcad.
MULPOL input only contains charge, dipole and quadrupole.

The program mpprop2mulfit converts a Molcas MpProp file to the MULPOL format (up to quadrupoles only). The conversion works well. Look in the code to see the required (strange) order of the dipoles and the Stones format

Results

Note that the results of this program are NOT idential to a standard ESP fit (owing to charge penetration is avoided). On the contrary, the two results are quite different.
Unfortunately, the results still depend strongly on the vdW radii and the integration limits.

My only tested example is for acetamid (CH3CONH2):

	MK	ChelpG	MulFit With bond centres			Mulfit Without bond centres			Original multipoles
Charges			Default	MK-like	ChelpG-like	Default	MK-like	ChelpG-like	Charges	c+dip	c+d+q
C1	-0.58	-0.43	-2.38	-2.68	-2.40	-1.05	-0.31	-0.97
C2	0.95	0.96	1.48	1.56	1.57	1.18	0.97	0.90
O	-0.65	-0.67	-0.74	-0.71	-0.74	-0.71	-0.70	-0.66
N	-1.13	-1.12	-1.16	-1.47	-1.48	-1.19	-1.16	-0.79
HC1	0.16	0.12	0.65	0.74	0.66	0.29	0.10	0.29
HC2	0.16	0.11	0.61	0.70	0.62	0.26	0.07	0.26
HC3	0.16	0.12	0.62	0.71	0.63	0.27	0.08	0.27
HN1	0.47	0.46	0.47	0.60	0.60	0.48	0.48	0.38
HN2	0.45	0.44	0.44	0.56	0.56	0.47	0.46	0.32
Dipole	0.017	0.014	0.014	0.021	0.013	0.054	0.065	0.066	0.185	0.001	0.001
Quadrupole	0.01	0.06	0.22	0.38	0.33	0.15	0.15	0.22	0.38	0.15	0.003
Octupole	1.8	2.2	2.8	3.4	2.0	1.8	2.1	1.6	0.4	0.4	0.4
Potential	0.0007	0.0009	0.0039	0.0046	0.0040	0.0013	0.0014	0.0020	0.0065	0.0016	0.0016

It can be seen that MulFit gives much worse results than the ESP charges (MK or CHELPG), especially for the q-pole and potential. However, the last three columns show that the original atomic multipoles do not reproduce the potential better than the MulFit. Thus, the difference is mainly caused by charge penetration.
The results with bond centres are obtained with a vdW radius of 1.6 Å.
The results are in Rhodopsin/Mulfit

Sample parameter file

With my suggested values for fitting of charges from MUPOL input:

TITLE
title
PRINT 2
VDW
6   1.75
CUMOVR 0
TOTMUL 3 1 0.0 0.0 0.0
REFRNK 2

The results are strongly sensitive to
INTLIM 1.2 3.0 - the integral thresholds (default values given). It affects only the integral limits, not the included atoms.
ADDRAD 0.0     - a constant added to the radii to define what atoms to include in the fit (default = 0 Å). It affects what atoms are included in the fit, not the integral limits, but it must be <= the lower INTLIM.

To avoid to fit any moment (even a charge) at a centre (e.g. at a bond centre), you can write (note that the second line is needed).
FITMOM
11 0
    0

Our testing indicates that a vdW radius of 1.6 Å gives best results for the bond centres
VDW
99   1.60

MK-like results are obtained with the following data
VDW
1   1.68
6   2.10
7   2.10
8   1.96
99   1.70
INTLIM 0.0 0.9

CHELP-G-like results are obtained with the following data
VDW
1   1.45
6   1.50
7   1.70
8   1.70
99   1.70
INTLIM 0.0 2.8

Sample MULPOL input file

The first line is needed (it marks that it is a MULPOL input)
The next 6 lines are ignored, but they need to be exactly 6.
Then, two lines per atom are read.
The first gives the atomic number and the xyz coordinates in A.
The second gives the multipole moments in the order: charge, dipole (11s, 10, 11c =y, z, x ), and quadrupole (22s, 21s, 20, 21c, 22c = xy, yz, zz, xz, xx-yy, Stone form)

DMASEM
THIS FILE CONTAINS DISTRIBUTED MULTIPOLES
CALCULATED FROM A SEMIEMPIRICAL WAVE-FUNCTION.
IT IS THE OUTPUT OF THE MULPOL PROGRAM AND
THE INPUT OF THE MULFIT PROGRAM
Two more lines are needed to read the file correctly
Two more lines are needed to read the file correctly
    6    0.0000    0.0000    0.0000
   -0.2416   -0.0098   -0.0017    0.0771    0.0088   -0.0086   -0.0429   -0.0013    0.1189
    6    1.5086    0.0000    0.0000
    0.2998   -0.0308   -0.0009    0.0984   -0.0721    0.0000    0.1273    0.0012   -0.0608
    8    2.1700    1.0585    0.0000
   -0.3748   -0.2286   -0.0001   -0.1394    0.2814   -0.0011    0.0063    0.0019   -0.1511
    7    2.1587   -1.2111   -0.0069
   -0.4426    0.0137   -0.0005   -0.0132   -0.0154    0.0037   -0.4528   -0.0030   -0.0175
    1   -0.3584    1.0474    0.1504
    0.1174    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1   -0.4020   -0.6409    0.8206
    0.0948    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1   -0.3800   -0.3781   -0.9797
    0.1002    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1    3.1476   -1.2342   -0.0113
    0.2288    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    1    1.6612   -2.0612   -0.0090
    0.2179    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

The Stone form is obtained in the following way (sqrt=sqrt(3)):
         result(1)=zz-0.5d0*(xx+yy)     ! 20
         result(2)=sqrt3*xz             ! 21c
         result(3)=sqrt3*yz             ! 21s
         result(4)=sqrt3*0.5d0*(xx-yy) ! 22c
         result(5)=sqrt3*xy             ! 22s

Sample CADPAC input file

Note that the first line is needed to start to read the input.
Next the program looks for 'x =' to find the coordinates and atom symbol
Then, it reads the lines with Q (but not |Q) are read.
Note that this is an old CADPAC 5.2 format file that is not correctly read by the original progam any longer, but it works after a few simple modifications of the program (see below).

                                             Distributed Multipole Analysis

Analysing the density matrix on section 7
(Default SCF section)

Mulliken Charge Densities
Basis Function Populations
     1     F s       1.99533
     2     F s       0.97089
     3     F x       0.91365
     4     F y       1.19279
     5     F z       1.19279
     6     F s       0.88213
     7     F x       0.63535
     8     F y       0.80558
     9     F z       0.80558
    10     F xx      0.06703
    11     F yy      0.02626
    12     F zz      0.02626
    13     F xy      0.00163
    14     F xz      0.00163
    15     F yz      0.00000
    16     H s       0.40820
    17     H s       0.07491

Atomic Charges
     1    F1          -0.51690
     2    H2           0.51690
Site 1    F1         x = 0.000000 y = 0.000000 z = 0.000000   Maximum
rank = 10   Relative radius = 1.000
                    Q00 = -0.510267
|Q1| =   0.132001 Q11c = -0.132001
|Q2| =   0.130714 Q20 = -0.065357 Q22c =   0.113202
|Q3| =   0.268791 Q31c =   0.164600 Q33c = -0.212498
|Q4| =   0.222455 Q40 = -0.083421 Q42c =   0.124356 Q44c = -0.164508
|Q5| =   0.177840 Q51c = -0.086097 Q53c =   0.092995 Q55c = -0.124766
|Q6| =   0.146851 Q60 =   0.045891 Q62c = -0.066502 Q64c =   0.072850
Q66c = -0.098639
|Q7| =   0.123988 Q71c =   0.051257 Q73c = -0.053267 Q75c =   0.058889
Q77c = -0.080253
|Q8| =   0.106588 Q80 = -0.029145 Q82c =   0.041802 Q84c = -0.043842
Q86c =   0.048783 Q88c = -0.066799
|Q9| =   0.092737 Q91c = -0.034021 Q93c =   0.034821 Q95c = -0.036757
Q97c =   0.041096 Q99c = -0.056481
|QA| =   0.081201 QA0 =   0.019983 QA2c = -0.028521 QA4c =   0.029381
QA6c = -0.031163 QA8c =   0.034970 QAAc = -0.048203

Site 2    H2         x = 1.732879 y = 0.000000 z = 0.000000   Maximum
rank = 10   Relative radius = 1.000
                    Q00 =   0.510267
|Q1| =   0.027821 Q11c =   0.027821
|Q2| =   0.043548 Q20 = -0.021774 Q22c =   0.037714
|Q3| =   0.063374 Q31c =   0.038809 Q33c = -0.050102
|Q4| =   0.061727 Q40 =   0.023148 Q42c = -0.034507 Q44c =   0.045648
|Q5| =   0.052207 Q51c = -0.025274 Q53c =   0.027300 Q55c = -0.036626
|Q6| =   0.041038 Q60 = -0.012824 Q62c =   0.018584 Q64c = -0.020358
Q66c =   0.027565
|Q7| =   0.030840 Q71c =   0.012749 Q73c = -0.013249 Q75c =   0.014648
Q77c = -0.019961
|Q8| =   0.022479 Q80 =   0.006147 Q82c = -0.008816 Q84c =   0.009246
Q86c = -0.010288 Q88c =   0.014088
|Q9| =   0.016028 Q91c = -0.005880 Q93c =   0.006018 Q95c = -0.006353
Q97c =   0.007103 Q99c = -0.009762
|QA| =   0.011239 QA0 = -0.002766 QA2c =   0.003948 QA4c = -0.004067
QA6c =   0.004313 QA8c = -0.004840 QAAc =   0.006672

Total multipoles referred to origin at x =   0.000000, y =    0.000000, z =

Compilation

Goto src and run compile.sh

You will get the following warning, which can be ignored:

rdat.f: In function `inpkey':
rdat.f:45: warning:
        ^       inpkey(cfl)+1
                1
rdat.f:379: (continued):
           function inpkey(fl)
                    2
Argument #1 (named `fl') of `inpkey' is one type at (2) but is some other type at (1) [info -f g77 M GLOBALS]

Tests

The formamid example gives a file for018 identical to formamid.sum, so I guess the program works properly. The formatmid.txt file cannot be compared to for017, because the former was removed by your mail server.

However, note that the results are wrong: You have to add another two title lines to make the program read all nine atoms. In the sample, only the last eight atoms are read and treated.

For the HF example, however, I get an error message:
> ln -fs hf.mlt for014
>mulfit
fmt: read unexpected character
apparent state: unit 14 named for014
last format: (15x,3(f10.6,5x))
lately reading sequential formatted external IO
Aborted (core dumped)

The end of the fort.30 file is:
Line:      2    H2           0.51690
Line,gen: Site 1    F1         x = 0.000000 y
NN = 1
Site 1    F1         x = 0.000000 y = 0.000000 z = 0.000000   M

To run this example (which has an old format of the CADPAC file, I had to do the following changes into the code (routine rcad.f):

14d13
<       logical cad5 ! UR 24/10-07
42d40
<       cad5=(index(line,'Site').ne.0) ! UR to identify type of input
102,106c100,101
<       if (cad5) then ! changed by UR
<          read (iref14,'(26x,3(f10.5,5x))') x,y,z    ! cadpac 5.2
<       else
<          read (iref14,'(15x,3(f10.6,5x))') x,y,z    ! cadpac 6.1
<       endif
---
> C     read (iref14,'(26x,3(f10.5,5x))') x,y,z    ! cadpac 5.2
>       read (iref14,'(15x,3(f10.6,5x))') x,y,z    ! cadpac 6.1