UT theoretical chemistry code forum

Dear all, I use 14 image to perform a neb running. After some inoic steps, the neb running is still not convergent. The following are the maximum forces on the 11th image:

NEB: Forces: max atom, RMS 0.631483 0.179357
NEB: Forces: max atom, RMS 0.440029 0.109492
NEB: Forces: max atom, RMS 0.398472 0.139530
NEB: Forces: max atom, RMS 0.286875 0.074952
NEB: Forces: max atom, RMS 0.385002 0.139893
NEB: Forces: max atom, RMS 0.294194 0.075594
NEB: Forces: max atom, RMS 0.304148 0.109699
NEB: Forces: max atom, RMS 0.448549 0.097984
NEB: Forces: max atom, RMS 0.550919 0.117665
NEB: Forces: max atom, RMS 0.344926 0.081073
NEB: Forces: max atom, RMS 0.282825 0.084797
NEB: Forces: max atom, RMS 0.182802 0.057848
NEB: Forces: max atom, RMS 0.252353 0.087619
NEB: Forces: max atom, RMS 0.148363 0.057552
NEB: Forces: max atom, RMS 0.252159 0.082422
NEB: Forces: max atom, RMS 0.156565 0.060299
NEB: Forces: max atom, RMS 0.161386 0.062582
NEB: Forces: max atom, RMS 0.109948 0.040971
NEB: Forces: max atom, RMS 0.105013 0.043827
NEB: Forces: max atom, RMS 0.089034 0.036491
NEB: Forces: max atom, RMS 0.089474 0.038830
NEB: Forces: max atom, RMS 0.122034 0.038181
NEB: Forces: max atom, RMS 0.126389 0.039811
NEB: Forces: max atom, RMS 0.096752 0.035446
NEB: Forces: max atom, RMS 0.109231 0.038393
NEB: Forces: max atom, RMS 0.103796 0.027472
NEB: Forces: max atom, RMS 0.150193 0.043565
NEB: Forces: max atom, RMS 0.147504 0.036362

You could find the maximum forces vibrate around the value of 0.1 and 0.2, and it does not descend to smaller value anymore.
So I want to know how to deal with this problem, whether I should take pains to wait more calculational iterations, or change the spring constant k (initially I set k=-5.0 in INCAR) , or other something?

Thanks a lot for your any advice.

Hi there,
I am new with NEB and am trying to find the MEP of extracting hydrogen from surface. I successfully relaxed structures for the initial and final configurations and once I ran "nebmake.pl" scripts it gave me the images in between, however when I ran vasp, in the OUTCAR file it complains that some of the atoms are too close to each other and thus it leads to the termnination of the job after some iterations.

Does anyone knows how to deal with the situation?
Any suggestion is appreciated,
Thanks in advance,
ali

@jzhao:

Try using a conservative optimizer to make sure you see systematic convergence. Using IBRION=3 and POTIM=0.1 is a very safe setting. It will not be the most efficient, but the force should systematically drop. I do not recommend changing the spring constant.

Also, check the energy. If the forces are fluctuating around 0.1 or 0.2 eV/Ang, this only indicates poor convergence if the energy is not changing. If the energy of the image is systematically changing, the calculation could be working properly.

Finally, we have found better convergence using our force-based conjugate gradient or FIRE optimizers. If the band is converging slowly, one of these optimizers could improve the convergence rate.

@aakbarz:

Ideally, we should add a function to push atoms apart in the nebmake.pl script. If there is enough demand, we'll add it (we use this in other codes).

To fix the problem with vasp, you can run a few very conservative optimization steps. Using IBRION=3 and POTIM=0.01 for 10 or 20 iterations is usually enough to bring the forces down to a reasonable 1-10 eV/Ang. At this point, POTIM can be raised to 0.1 or a more aggressive optimizer can be used. If you are using our new optimizers, you can also constrain the optimization by choosing a small value of MaxMove.

Thanks in advance for your feedback.
Just another question, what kind of function should one add to push the atoms apart in the nebmake.pl ?

In which code do you use the function?

thanks again for the help,
best,

Hello aakbarz, the function we used is part of the akmc script. I just wrapped it into a simple perl script so that you can use it independently.
uasge:
1. copy the script lines into a file and name the file as "neb_pushapart.pl"
2. put the file into the same directory as "nebmake.pl" resides in.
3. make the file an executable command (for example, chmod 750 neb_pushapart.pl)
4. type "rehash";
5. go into the neb directory
6. type "pushapart.pl 0.3" and the allowed distance between any two atoms in the system is no more than 0.3 Angstrom; otherwise, a spring force would be applied to those involved atoms and push them apart little by little. You can change the value as you like, but it should not be a big value; 0.1-0.3 should be enough.

if you simply type "pushapart.pl", the default allowed distance between two atoms is no more than 0.1 Angstrom.

The script will go into each image folder and save the old POSCAR as POSCAR_orig. The POSCAR file will be rewritten.

Please let me know if it works out. OK, the followings are the script lines.

eval '(exit $?0)' && eval 'exec perl -S $0 ${1+"$@"}' && eval 'exec perl -S $0 $argv:q' if 0;
#;-*- Perl -*-
use FindBin qw($Bin);
use lib "$Bin";
use Vasp;

$rcut=0.1;
if(@ARGV>0) {$rcut=$ARGV[0];}
@folders=`ls -d [0-9][0-9]`;
#print "@folders\n";
$num=@folders-1;
$POSCAR="POSCAR";
for($i=1;$i<$num;$i++){
chomp($targetfolder=$folders[$i]);
$targetfolder=~s/^\s+//;
print "working on $targetfolder...\n";
($coordinates,$basis,$lattice,$num_atom,$total_atoms,$selectiveflag,$selective,$description)=read_poscar($POSCAR);
spring_relaxation($coordinates,$basis,$lattice,$total_atoms,$selective,$rcut);
system "cd $targetfolder; mv $POSCAR POSCAR_orig";
$outputfilename=$targetfolder."/POSCAR";
write_poscar($coordinates,$basis,$lattice,$num_atom,$total_atoms,$selectiveflag,$selective,$description,$outputfilename);
}
sub spring_relaxation{
my ($R,$basis,$lattice,$totalatoms,$selective,$cutoff)=@_;
my $epsilon=0.1;
my $stepmax=10000;
my $drmax=0.1;
my @line=();
my ($i,$j,$Rijx,$Rijy,$Rijz,$Ax,$Ay,$Az,$Ax2,$Ay2,$Az2,$Ax2_,$Ay2_,$Az2_,$Rij);
my ($force,$vector1,$vector2,$vector3,$fij,$phiij,$step,$goodsound,$converged,$frozen,$dr);
# direct to cartesian
$R=dirkar($R,$basis,$lattice,$totalatoms);
# figure out box dimension
for($j=0;$j<3;$j++) {
$vector1->[0][$j]=$basis->[$j][0];
$vector2->[0][$j]=$basis->[$j][1];
$vector3->[0][$j]=$basis->[$j][2];
}
$Ax=magnitude($vector1,1);
$Ay=magnitude($vector2,1);
$Az=magnitude($vector3,1);
$Ax2=0.5*$Ax;
$Ay2=0.5*$Ay;
$Az2=0.5*$Az;
$Ax2_=-0.5*$Ax;
$Ay2_=-0.5*$Ay;
$Az2_=-0.5*$Az;
$Ax_=-1.0*$Ax;
$Ay_=-1.0*$Ay;
$Az_=-1.0*$Az;
# figure out frozen flags for each degree of freedom
for($i=0;$i<$totalatoms;$i++) {
@line=split(/\s+/,$selective->[$i]);
for($j=0;$j<3;$j++) {
$frozen->[$i][$j]=$line[$j];
}
}
# start the steepest descent loop
$converged=0;
for($step=0;$step<$stepmax;$step++) {
if($converged) {last;}
#print "Step No. $step\n";
for($i=0;$i<$totalatoms;$i++) {
for($j=0;$j<3;$j++) {
$force->[$i][$j]=0.0;
}
}
$goodsound=1;
for($i=0;$i<$totalatoms-1;$i++) {
for($j=$i+1;$j<$totalatoms;$j++) {
$Rijx=$R->[$i][0] - $R->[$j][0];
$Rijy=$R->[$i][1] - $R->[$j][1];
$Rijz=$R->[$i][2] - $R->[$j][2];
if($Rijx < $Ax2_){$Rijx+=$Ax;}
if($Rijy < $Ay2_){$Rijy+=$Ay;}
if($Rijz < $Az2_){$Rijz+=$Az;}
if($Rijx > $Ax2){$Rijx-=$Ax;}
if($Rijy > $Ay2){$Rijy-=$Ay;}
if($Rijz > $Az2){$Rijz-=$Az;}
$Rij=sqrt($Rijx*$Rijx+$Rijy*$Rijy+$Rijz*$Rijz);
if($Rij < $cutoff){
$goodsound=0;
($fij,$phiij)=linear_repulsion_pot($Rij,$cutoff);
#($fij,$phiij)=exp_repulsion_pot($Rij,$cutoff);
$force->[$i][0]=$force->[$i][0]-$fij*$Rijx/$Rij;
$force->[$i][1]=$force->[$i][1]-$fij*$Rijy/$Rij;
$force->[$i][2]=$force->[$i][2]-$fij*$Rijz/$Rij;
$force->[$j][0]=$force->[$j][0]+$fij*$Rijx/$Rij;
$force->[$j][1]=$force->[$j][1]+$fij*$Rijy/$Rij;
$force->[$j][2]=$force->[$j][2]+$fij*$Rijz/$Rij;
}
}
}
if($goodsound) {
$converged=1;
next;
}
for($i=0;$i<$totalatoms;$i++) {
for($j=0;$j<3;$j++) {
#if($frozen->[$i][$j] == 1) {next;}
if(lc($frozen->[$i][$j]) eq "f") {next;}
$dr=$epsilon*$force->[$i][$j];
if(abs($dr) > $drmax) {$dr=$drmax*$dr/abs($dr);}
#print "$i: $j -- dr: $dr\n";
$R->[$i][$j]+=$dr;
}
while ($R->[$i][0] < 0) {$R->[$i][0]+=$Ax;}
while ($R->[$i][1] < 0) {$R->[$i][1]+=$Ay;}
while ($R->[$i][2] < 0) {$R->[$i][2]+=$Az;}
while ($R->[$i][0] > $Ax) {$R->[$i][0]-=$Ax;}
while ($R->[$i][1] > $Ay) {$R->[$i][1]-=$Ay;}
while ($R->[$i][2] > $Az) {$R->[$i][2]-=$Az;}
}
}
if($converged == 1) {
print "fully relaxed after $step (1: trivial) steps\n";
}else{
print "Not fully relaxed after $stepmax steps, but we can't wait, so move on.\n";
}
# Cartesian to direct
$R=kardir($R,$basis,$lattice,$totalatoms);
}

sub linear_repulsion_pot{ # E=kx-b, k< 0
my ($r,$cutoff)=@_;
my $slope=-0.1;
my ($energy, $fij);
my $intercept=$slope*$cutoff;
$energy=$slope*$r-$intercept;
$fij=$slope;
return ($fij,$energy);
}

Is this a hydrogen molecule that you removed from the surface? Then make sure that the two different H atoms do no cross in the NEB, i.e. they are in the same order in the two original POSCAR files. It is also easy to use the neb_movie.pl script to generate a movie (.xyz file) of the initial linear NEB. Any crossing should be obvious from there.

Hi there,
Thanks a lot for the script. I tried your tips step by step but when I ran "neb_pushapart.pl", it does not give me e.g. POSCAR.orig, thus I am not sure if it is working or not! and also it gives me the following masseg:

working on 01...
fully relaxed after 1 (1: trivial) steps
Illegal division by zero at /Users/vtstscripts/Vasp.pm line 447.

Is anything wrong with my POSCAR's?
could you please check the script, I can also post my POSCAR's files for you to check.

thanks in advance,
best,
ali

Could you post your POSCAR files in the neb folder 01 and 02?

[quote="xulijun2"]Could you post your POSCAR files in the neb folder 01 and 02?[/quote]


Here is POSCAR in 01

libh4
1.00000000000000
7.2169999999999996 0.0000000000000000 0.0000000000000000
0.0000000000000000 4.3893000000000004 0.0000000000000000
0.0000000000000000 0.0000000000000000 26.5072000000000010
8 32 8
Selective dynamics
Direct
0.2100382923053770 0.2581497752031720 0.0923854573080831 T T T
0.2984086249417150 0.2524272274114450 0.3118401708422030 T T T
0.1877054039749030 0.7507577441460070 0.2333673032042970 T T T
0.0652961694267873 0.7354698065533470 0.4111132981421550 T T T
0.7213155605011040 0.8372807314657250 0.1653699965540480 T T T
0.6164402825840870 0.8520495082508009 0.3848134437981300 T T T
0.8879889010775560 0.1630993656214910 0.1133075982269390 T T T
0.7815084253766860 0.3448640251666790 0.2944079217622640 T T T
0.8099690078645310 0.2442332166206430 0.2133440858752850 T T T
0.0091340000182285 0.9566225613069220 0.4422429331252220 T T T
0.5841957536632170 0.7012185110451060 0.1060160889815340 T T T
0.6589286920332280 0.7079466670006130 0.2918477347560890 T T T
0.0496707918501826 0.1353124960658720 0.0847153487101622 T T T
0.0754893499526972 0.7140080788654610 0.3151798260471140 T T T
0.4404348655065320 0.3417457259890780 0.1897317446524620 T T T
0.3430552676326000 0.3884196014342930 0.3714082038547490 T T T
0.4846893637785000 0.2059537597989090 0.0673560195084448 T T T
0.3792026265530740 0.1637107680949940 0.2746779060064670 T T T
0.0945183760053254 0.7482123997070720 0.1946214809089350 T T T
0.2075684840754590 0.6627217469704620 0.3853872802891090 T T T
0.6765624809255540 0.9273759381770220 0.2023311859976500 T T T
0.5639371288224591 0.8683329028825570 0.4083605955515210 T T T
0.9082227640003500 0.2832563774807260 0.1476940342544920 T T T
0.8346193772474140 0.3492811568565910 0.3327919060645070 T T T
0.1411854378048790 0.0648485739851452 0.0863028972864299 T T T
0.1825566512918900 0.1532193207308130 0.3239133604310570 T T T
0.2906977092379820 0.5300281507676060 0.2332427833491690 T T T
0.1111242564848850 0.6160043187702510 0.4343082958437140 T T T
0.2207058886727680 0.5665712066511400 0.1331891493595240 T T T
0.1328063646091380 0.5010224104258040 0.3101272403793590 T T T
0.3650849181166010 0.0645396305463990 0.1904107027038360 T T T
0.2843217711641940 0.9118353598908731 0.3999384827145960 T T T
0.7961367711065320 0.9353488243416170 0.1532672196490310 T T T
0.8805572413015030 0.0694205340098151 0.3826778863287070 T T T
0.8892869978276470 0.5437798220407610 0.9300623209970670 T T T
0.7179851199990140 0.5665216234072050 0.2945857046517730 T T T
0.8743734498076420 0.5526557817098380 0.0596530680539828 T T T
0.7868766174904490 0.5255394909128980 0.3956906589636300 T T T
0.6939228421345400 0.1769013339219610 0.0772422452353183 T T T
0.6113088746641320 0.0179811544799975 0.2854924085472990 T T T
0.8286966671182870 0.8690730590958680 0.0546885340054715 T T T
0.1452053976335960 0.2532350484031340 0.2317927467211690 T T T
0.3362561360164660 0.7398534010618580 0.1523944033414320 T T T
0.4569835794399070 0.7550097645432100 0.3661655111956460 T T T
0.8046103500687640 0.8406470871760310 0.2486353584508340 T T T
0.1289123326004710 0.2285318141409080 0.4214733980867180 T T T
0.5700383290226581 0.3368516293489560 0.1613200008117420 T T T
0.6876855774686210 0.1578174560031020 0.3747168755225730 T T T


POSCAR in 02

libh4
1.00000000000000
7.2169999999999996 0.0000000000000000 0.0000000000000000
0.0000000000000000 4.3893000000000004 0.0000000000000000
0.0000000000000000 0.0000000000000000 26.5072000000000010
8 32 8
Selective dynamics
Direct
0.1181728742615750 0.2744331295429990 0.0765003780196594 T T T
0.3015080846855160 0.2450998428795760 0.2617638462291650 T T T
0.1861907811897580 0.7464954648761970 0.2341252782168300 T T T
0.9693257270618220 0.7357895636341421 0.3461063690965030 T T T
0.7432510528728000 0.9362004693733610 0.1931691605759980 T T T
0.5407801312684200 0.9476554222312930 0.3785858513144830 T T T
0.9443834438622790 0.0632625216713336 0.2034945860504960 T T T
0.7557419620832631 0.4421083812249900 0.3217883407482990 T T T
0.7117573795362450 0.2386850762763260 0.1968139448908700 T T T
0.9846310008290190 0.0259217797313198 0.3915938956317930 T T T
0.5830985871304220 0.7071980089527660 0.1080313408187190 T T T
0.7269645750677900 0.6656585976838800 0.2302927541565080 T T T
0.1396503066430630 0.1590250462539800 0.1630365399185790 T T T
0.0628258851201409 0.6763455532357060 0.3606313512685660 T T T
0.4778614396847730 0.4374444587223110 0.2335196624605800 T T T
0.2779650642460400 0.4740679804996330 0.3470256783800640 T T T
0.5775035062898159 0.1831132399759670 0.0654736178474593 T T T
0.3829528633290520 0.1793812964806420 0.2237684300901180 T T T
0.0927254646790203 0.7454604962396230 0.1954130881030810 T T T
0.3029253709139470 0.6810814021373230 0.3365707882998880 T T T
0.7327306368006320 0.0037098051593674 0.2306514882931620 T T T
0.5255247037094650 0.9619764624463660 0.3866703738847550 T T T
0.9360239208220240 0.2122245422125320 0.2303705234836500 T T T
0.7691450100483480 0.4457647453298180 0.3597772431234740 T T T
0.0850274116791354 0.1060701229587480 0.0623109947866141 T T T
0.1887785084070120 0.2254613591773730 0.2712616259330290 T T T
0.2885176992879150 0.5252665518613210 0.2338092154804630 T T T
0.0271560435928677 0.7034231703760910 0.3632556155633680 T T T
0.2362830343557220 0.6636621955504650 0.1590043921546500 T T T
0.0462600680340728 0.5057895322983870 0.2706189469483180 T T T
0.4440113105304060 0.1463818366616250 0.1483988139014020 T T T
0.2658714966242850 0.9337446621461610 0.3280664497464580 T T T
0.7741853160443950 0.9542613126766850 0.1840335146586720 T T T
0.9633675737306360 0.1658020870815090 0.3769986487517800 T T T
0.9812850854437440 0.6184980450682500 0.8662888119216490 T T T
0.7322925831375220 0.6651236393198370 0.3223878294768310 T T T
0.9735645878726020 0.6036200103812061 0.9694909463895630 T T T
0.7859217177419920 0.5233270697466070 0.3988028651453420 T T T
0.6974421737167090 0.2463324541873320 0.1270471574983740 T T T
0.5123872415561550 0.0159592714069561 0.3035869385427010 T T T
0.7861612584304321 0.9677052801046170 0.0660859750534826 T T T
0.1368857593580660 0.2447583864733500 0.1861855080417740 T T T
0.3376051737475050 0.7354259639173220 0.1544218714118910 T T T
0.5518560810487190 0.7493824324241700 0.3319395924939390 T T T
0.7666770900477360 0.9405544097648429 0.2751314903933770 T T T
0.1356935276070890 0.2296409924510170 0.3871816409325390 T T T
0.5081896906469440 0.4363637900717310 0.2235772854726430 T T T
0.7158147353476350 0.0610674896490764 0.4004683925467130 T T T

thanks a lot,

Hey there, please try the following script. I hope it will work this time.

eval '(exit $?0)' && eval 'exec perl -S $0 ${1+"$@"}' && eval 'exec perl -S $0 $argv:q' if 0;
#;-*- Perl -*-
use FindBin qw($Bin);
use lib "$Bin";
use Vasp;

$rcut=0.1;
if(@ARGV>0) {$rcut=$ARGV[0];}
chomp($folder=`ls -d [0-9][0-9]`);
$folder=~s/^\s+//;
@folders=split(/\s+/,$folder);
$num=@folders-1;
for($i=1;$i<$num;$i++){
$targetfolder=$folders[$i];
$dummy=substr($targetfolder, -1 , 1);
if($dummy eq "\/") {chop($targetfolder);}
$POSCAR=$targetfolder."/POSCAR";
if(!(-e $POSCAR)) {die "$POSCAR does not exist.\n";}
print "working on $targetfolder...\n";
($coordinates,$basis,$lattice,$num_atom,$total_atoms,$selectiveflag,$selective,$description)=read_poscar($POSCAR);
spring_relaxation($coordinates,$basis,$lattice,$total_atoms,$selective,$rcut);
system "cd $targetfolder; mv POSCAR POSCAR_orig";
$outputfilename=$targetfolder."/POSCAR";
write_poscar($coordinates,$basis,$lattice,$num_atom,$total_atoms,$selectiveflag,$selective,$description,$outputfilename);
}
sub spring_relaxation{
my ($R,$basis,$lattice,$totalatoms,$selective,$cutoff)=@_;
my $epsilon=0.1;
my $stepmax=10000;
my $drmax=0.1;
my @line=();
my ($i,$j,$Rijx,$Rijy,$Rijz,$Ax,$Ay,$Az,$Ax2,$Ay2,$Az2,$Ax2_,$Ay2_,$Az2_,$Rij);
my ($force,$vector1,$vector2,$vector3,$fij,$phiij,$step,$goodsound,$converged,$frozen,$dr);
# direct to cartesian
$R=dirkar($R,$basis,$lattice,$totalatoms);
# figure out box dimension
for($j=0;$j<3;$j++) {
$vector1->[0][$j]=$basis->[$j][0];
$vector2->[0][$j]=$basis->[$j][1];
$vector3->[0][$j]=$basis->[$j][2];
}
$Ax=magnitude($vector1,1);
$Ay=magnitude($vector2,1);
$Az=magnitude($vector3,1);
$Ax2=0.5*$Ax;
$Ay2=0.5*$Ay;
$Az2=0.5*$Az;
$Ax2_=-0.5*$Ax;
$Ay2_=-0.5*$Ay;
$Az2_=-0.5*$Az;
$Ax_=-1.0*$Ax;
$Ay_=-1.0*$Ay;
$Az_=-1.0*$Az;
# figure out frozen flags for each degree of freedom
for($i=0;$i<$totalatoms;$i++) {
@line=split(/\s+/,$selective->[$i]);
for($j=0;$j<3;$j++) {
$frozen->[$i][$j]=$line[$j];
}
}
# start the steepest descent loop
$converged=0;
for($step=0;$step<$stepmax;$step++) {
if($converged) {last;}
#print "Step No. $step\n";
for($i=0;$i<$totalatoms;$i++) {
for($j=0;$j<3;$j++) {
$force->[$i][$j]=0.0;
}
}
$goodsound=1;
for($i=0;$i<$totalatoms-1;$i++) {
for($j=$i+1;$j<$totalatoms;$j++) {
$Rijx=$R->[$i][0] - $R->[$j][0];
$Rijy=$R->[$i][1] - $R->[$j][1];
$Rijz=$R->[$i][2] - $R->[$j][2];
if($Rijx < $Ax2_){$Rijx+=$Ax;}
if($Rijy < $Ay2_){$Rijy+=$Ay;}
if($Rijz < $Az2_){$Rijz+=$Az;}
if($Rijx > $Ax2){$Rijx-=$Ax;}
if($Rijy > $Ay2){$Rijy-=$Ay;}
if($Rijz > $Az2){$Rijz-=$Az;}
$Rij=sqrt($Rijx*$Rijx+$Rijy*$Rijy+$Rijz*$Rijz);
if($Rij < $cutoff){
$goodsound=0;
($fij,$phiij)=linear_repulsion_pot($Rij,$cutoff);
#($fij,$phiij)=exp_repulsion_pot($Rij,$cutoff);
$force->[$i][0]=$force->[$i][0]-$fij*$Rijx/$Rij;
$force->[$i][1]=$force->[$i][1]-$fij*$Rijy/$Rij;
$force->[$i][2]=$force->[$i][2]-$fij*$Rijz/$Rij;
$force->[$j][0]=$force->[$j][0]+$fij*$Rijx/$Rij;
$force->[$j][1]=$force->[$j][1]+$fij*$Rijy/$Rij;
$force->[$j][2]=$force->[$j][2]+$fij*$Rijz/$Rij;
}
}
}
if($goodsound) {
$converged=1;
next;
}
for($i=0;$i<$totalatoms;$i++) {
for($j=0;$j<3;$j++) {
#if($frozen->[$i][$j] == 1) {next;}
if(lc($frozen->[$i][$j]) eq "f") {next;}
$dr=$epsilon*$force->[$i][$j];
if(abs($dr) > $drmax) {$dr=$drmax*$dr/abs($dr);}
#print "$i: $j -- dr: $dr\n";
$R->[$i][$j]+=$dr;
}
while ($R->[$i][0] < 0) {$R->[$i][0]+=$Ax;}
while ($R->[$i][1] < 0) {$R->[$i][1]+=$Ay;}
while ($R->[$i][2] < 0) {$R->[$i][2]+=$Az;}
while ($R->[$i][0] > $Ax) {$R->[$i][0]-=$Ax;}
while ($R->[$i][1] > $Ay) {$R->[$i][1]-=$Ay;}
while ($R->[$i][2] > $Az) {$R->[$i][2]-=$Az;}
}
}
if($converged == 1) {
print "fully relaxed after $step (1: trivial) steps\n";
}else{
print "Not fully relaxed after $stepmax steps, but we can't wait, so move on.\n";
}
# Cartesian to direct
$R=kardir($R,$basis,$lattice,$totalatoms);
}

sub linear_repulsion_pot{ # E=kx-b, k< 0
my ($r,$cutoff)=@_;
my $slope=-0.1;
my ($energy, $fij);
my $intercept=$slope*$cutoff;
$energy=$slope*$r-$intercept;
$fij=$slope;
return ($fij,$energy);
}

Hi there,
thanks a lot for your help, this one works.