[Molpro-user] Calculation of Potential III - SOLVED

Fri Sep 12 10:38:38 BST 2008

Hello All!

I have found a workaround for the problem and since similar thresholds  
may occur in other parts of Molpro as well I'd just like to post the  
solution. Maybe it's helpful for somebody stumbling over the question  
in the archive later on.

Two limits apply to PROPERTY:
- Per invocation in an execution block only 30 points can be  
calculated ("Too many operators" in propa.f thereafter)
- For each job at most 148 points can be calculated in all combined  
PROPERTY invocations ("Too many operators" in read_op.f afterwards)

The trick is to calculate the wavefunction and integrals, save the  
restart data (.wfu, .int & .aux) and calculate all needed points in  
portions of at most 148 using always a *fresh* copy of the scratch  
files. That is, copy the three files to new names, run a "restarted"  
job with e.g. 5 PROPERTY blocks with 29 points each, then delete the  
scratch files and create a new copy for the next 5 * 29 points. When  
restarting the job the files have to be in the same directory, copying  
the scratch files to a sub folder is not enough ("Recomputing  
integrals since basis changed" and "Illegal geometry input - no atoms").

Example input follows below. Hope this helps anyone!

Best regards,

Ben

=================================================
Example Input for the Property Calculations
=================================================

***, ELECPOT
file, 1, scratch.copy.int
file, 2, scratch.copy.wfu
file, 3, scratch.copy.aux

{matrop
read, savmat, DEN, charge, file=density.mat
save, savmat, 2140.2, DENSITY}

{property
density,2140.2
pot,,  -5.66918,  -7.55890,  -4.72431
pot,,  -5.66918,  -7.55890,  -3.77945
pot,,  -5.66918,  -7.55890,  -2.83459
  [ etc., etc., 23 other points here ]
pot,,  -5.66918,  -5.66918,   5.66918
pot,,  -5.66918,  -5.66918,   6.61404
pot,,  -5.66918,  -4.72431,  -6.61404}

[etc, 5 blocks like the above]

=================================================

On 10 Sep 2008, at 10:52, Benjamin Bulheller wrote:

> Hello All!
>
> Regarding my problem below I have discovered a new threshold in  
> Molpro which puzzles me.
>
> I can only calculate the electrostatic potential at 30 points per  
> calculation before it hits the maximum value 30 of the MXOPER  
> variable in propa.f. I then determined and saved the integrals in  
> one job and ran individual jobs re-using these data with an  
> evaluation of 30 points per run.
>
> This works for 148 points (as I found before) and then it hits the  
> maximum value of the variable MAXOPER (note the additional A, it is  
> actually set to 150 what is reached during a loop). This happens in  
> readop.f (routine write_op) and Molpro dies with the same error as  
> before, "Too many operators".
>
> Thanks a lot for any ideas or hints!
>
> Cheers,
>
> Ben
>
>
>
>
> Begin forwarded message:
>
>> From: Benjamin Bulheller <pcxbb2 at nottingham.ac.uk>
>> Date: 3 September 2008 17:44:06 BST
>> To: molpro-user at molpro.net
>> Subject: Calculation of Potential
>>
>> Hi!
>>
>> I need to calculate the electrostatic potential of a compound at  
>> several thousand points around the molecule. This has been done  
>> before (with Molpro 2002.6) using the following code:
>>
>>   property
>>   density,2100.2
>>   pot,,-12.28322, -2.83459, -0.94486
>>   pot,,-12.24522, -0.94486,  1.88973
>>   pot,,-12.67322,  1.88973,  0.94486
>>   pot,,-11.33836, -2.83459, -3.77945
>>   pot,,-10.12736,  2.83459,  2.83459
>>   etc...
>>
>> In 2002.6 this worked with up to 2500 "pot" statements while in  
>> 2006.1 only 30 points are possible and I get a "Too many operators"  
>> error with more than that. This happens in propa.f caused by the  
>> MXOPER variable.
>>
>> Using a slightly crude workaround I can wrap the pot statement in  
>> nested do loops and for some reason Molpro will then compute 148  
>> points and fail thereafter.
>>
>> Am I doing something wrong there in general? What puzzles me is  
>> that it obviously must have worked in 2002.6 but stopped working in  
>> 2006.1 (without changes in the concerned code section), is there  
>> any explanation for that?
>>
>> Thank you very much in advance,
>>
>> Ben
>>
>>
>>
>> =================================================
>> Example Input for H2
>> =================================================
>>
>> print,basis
>> gprint,orbitals
>> geomtyp=xyz
>>
>> geometry={
>> nosym
>> 2
>> Title
>> H 0.0 0.0 0.0
>> H 0.0 0.0 1.0
>> }
>>
>> basis={
>> ! HYDROGEN   [7s,3p] -> (2s,1p)  ANO-S
>> s,H , 
>> 82.636374,12.409558,2.8238540,0.79767000,0.25805300,0.08989100,0.03146200
>> c,1.7, 0.0025641, 0.0192126, 0.0950512, 0.3018678, 0.4874459,  
>> 0.23547943, 0.01577776
>> c,1.7,-0.0046008,-0.0324447,-0.1710737,-0.5048672,-0.0948920,  
>> 0.37085551, 0.65565134
>> p,H ,1.6625000,0.41560000,0.10390000
>> c,1.3, 0.3012771, 0.6800574, 0.1985627
>> }
>>
>> int;
>> hf;
>> property
>> density,2100.2
>> dm,,0.0,0.0,0.0;
>>
>> property
>> density,2100.2
>> pot,,-12.28322, -2.83459, -0.94486
>> pot,,-12.24522, -0.94486,  1.88973
>> pot,,-12.67322,  1.88973,  0.94486
>> pot,,-11.33836, -2.83459, -3.77945
>> pot,,-10.12736,  2.83459,  2.83459
>> pot,, -9.33836,  2.83459,  3.77945
>> pot,, -8.25736,  3.77945,  1.88973
>> ! plus 1000 and more additional lines for other points
>>
>>
>> =================================================
>
>
>
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.molpro.net/pipermail/molpro-user/attachments/20080912/1a040f1f/attachment.html>