#!/bin/sh

##########################################################################
# A molecular dynamics simulation of a sodium ion in a small water cluster.
# The oxygen atoms of the water molecules are constrained to be within 0.45 nm of the central ion.
#
# The constraint is a harmonic function U(r)=kappa*(r-0.45)**2.
# The constrained value r is the maximum distance out of the set of Na-O distances.
# It is a continuous function that approximates the maximum out of a set of values.
# The parameter beta determines the gradient of the function.
# Choosing a small beta value results in a better match between the actual
# and approximate maximum values, but also in a less continuous function.
##########################################################################

$AMSBIN/ams << eor

Task MolecularDynamics

MolecularDynamics
    NSteps 1000
    TimeStep 0.5
    InitialVelocities
        Temperature 300
    End
    Thermostat
        Type Berendsen
        Temperature 300
        Tau 100
    End
    Trajectory
        SamplingFreq 1
    End
    Plumed
        Input 
GROUP ATOMS=2,5,8,11,14,17,20,23,26 LABEL=oxygens
DISTANCES GROUPA=1 GROUPB=oxygens MAX={BETA=0.05} LABEL=wallcv
UPPER_WALLS ARG=wallcv.max AT=0.45 KAPPA=5000.0
PRINT ARG=wallcv.max STRIDE=1 FILE=COLVAR
        End
    End
End

system
  Atoms
             Na       0.0000000000       0.0000000000       0.0000000000
              O       1.7119252402      -2.0080034855       1.0155335456
              H       2.5731328659      -2.4678020432       1.2672765294
              H       1.4480244399      -2.5303870208       0.1729165635
              O       0.6300542373       0.5761441595       2.8226448718
              H       1.2391579321      -0.0089304214       2.3349735783
              H       0.9338206205       1.4729694268       2.4878846070
              O       0.1647045543       2.0208702482      -2.3249601218
              H      -0.7307215954       2.0625495942      -2.7797813945
              H       0.2097075116       1.0635963204      -2.1125352562
              O      -1.7025066874      -1.0780652110      -3.0700032660
              H      -2.3616744611      -0.4381207399      -3.4654803274
              H      -1.8354437418      -1.8757732529      -3.6203721171
              O       3.0658874271       0.7135749531      -2.1674092521
              H       3.9052706085       1.2362759268      -2.1017219618
              H       2.5231505825       1.3717470880      -2.6845066156
              O      -2.8805081861       2.6302885978      -0.3532568156
              H      -3.8762635503       2.6712519969      -0.4669367452
              H      -2.6280204633       2.2152595735      -1.2369080650
              O      -3.5076002124      -1.7860735453      -0.4319799526
              H      -3.8213733752      -1.9705946207      -1.3583130388
              H      -2.8954116001      -1.0294541930      -0.5563055520
              O       2.9868579456       2.5941122968       0.9208652453
              H       2.2983746780       3.1353478114       0.4089961955
              H       3.1905202107       1.8321429488       0.2760998864
              O      -2.9007338617       1.4515370844       2.5612078673
              H      -2.1648641984       1.4189029227       1.8879229015
              H      -3.6017400767       1.9414536506       2.1086226997
  End
  BondOrders
     20 21 1.0
     20 22 1.0
     26 27 1.0
     26 28 1.0
     23 24 1.0
     23 25 1.0
     2 3 1.0
     2 4 1.0
     11 12 1.0
     11 13 1.0
     17 18 1.0
     17 19 1.0
     14 15 1.0
     14 16 1.0
     5 6 1.0
     5 7 1.0
     8 9 1.0
     8 10 1.0
  End
  Charge 1.0
End

Engine ForceField
  Type UFF
EndEngine

RNGSeed 1

eor

##########################################################################
# The constrained values are stored in the file named COLVAR.
# Below we use Python to read the COLVAR file, and print the average, maximum, and minimum value.
##########################################################################

$AMSBIN/amspython << eor
infile = open('COLVAR')
lines = infile.readlines()
infile.close()

values = [float(line.split()[1]) for line in lines[1:]]
print()
print('Average value constraint: ',sum(values)/len(values))
print('Minimum value constraint: ',min(values))
print('Maximum value constraint: ',max(values))
eor