Example: Measuring the viscosity using NEMD

Download MD_shear_viscosity.run

#!/bin/sh

# This example performs a NEMD simulation of a water box under a shear acceleration profile.
# The acceleration profile is cosine-shaped and oriented along the Z axis (so that atoms with Z=0 are
# pulled in one direction, atoms in the middle of the box are pulled the other way). The flow direction
# is aligned with the X axis.
#
# This periodic perturbation induces a similarly cosine-shaped flow velocity profile in the liquid.
# The steady-state average amplitude of the flow velocity is then used to calculate the viscosity.

"$AMSBIN/ams" << eor
# Set the RNG seed to make results reproducible for regression testing, delete this line for real simulations
RNGSeed 836841470 1678420146 -1718043284 -369569254 -2053894812 -1628712594 -240676301 -356441007

Task MolecularDynamics
MolecularDynamics
    NSteps 24000
    TimeStep 1
    InitialVelocities
        Type Random
        Temperature 298.15
    End
    Thermostat
        Type NHC
        Temperature 298.15
        Tau 1000
    End
    CosineShear
        Enabled Yes
        Acceleration 0.000004522108063030400
    End
End
System
    Atoms
        O 0.8244444545793177 4.544838021389278 7.080687151044827 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.3547827848924093 4.719510963847915 7.916536224904227 ForceField.Type=HW ForceField.Charge=0.417
        H 1.355721739229628 5.344893665481732 6.975525637919818 ForceField.Type=HW ForceField.Charge=0.417
        O 4.824211318207558 0.8069313216022225 7.759825120988862 ForceField.Type=OW ForceField.Charge=-0.834
        H 4.514558713823999 1.280197293066311 8.558923240175481 ForceField.Type=HW ForceField.Charge=0.417
        H 4.033391063195659 0.7201928228115005 7.187127164561153 ForceField.Type=HW ForceField.Charge=0.417
        O 2.934572969015879 0.626831742380179 5.865967452739235 ForceField.Type=OW ForceField.Charge=-0.834
        H 2.151893964459041 1.137329352693095 6.027000900498897 ForceField.Type=HW ForceField.Charge=0.417
        H 3.263433084806719 1.048186940192297 5.040531567587069 ForceField.Type=HW ForceField.Charge=0.417
        O 1.206430947463354 1.856002602957112 9.003428097539595 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.7495045796559563 2.685732546074114 9.186271372016726 ForceField.Type=HW ForceField.Charge=0.417
        H 0.5732201719526651 1.426252646404271 8.428676649867203 ForceField.Type=HW ForceField.Charge=0.417
        O 4.00858752820043 4.029996794081304 5.583951652401872 ForceField.Type=OW ForceField.Charge=-0.834
        H 4.6767739512274 4.347838022636921 6.162932810282044 ForceField.Type=HW ForceField.Charge=0.417
        H 3.285915055183219 4.655218377978755 5.758534605518392 ForceField.Type=HW ForceField.Charge=0.417
        O 3.563629518129562 1.371347494469859 12.8855840322333 ForceField.Type=OW ForceField.Charge=-0.834
        H 3.213171669104573 0.4819670684498465 12.77900015824002 ForceField.Type=HW ForceField.Charge=0.417
        H 3.044964414139814 1.734726370517266 13.59851543041306 ForceField.Type=HW ForceField.Charge=0.417
        O 0.3789251709796683 1.127496658856301 13.2115287736376 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.6813483419435356 1.239300575078584 14.13801389632943 ForceField.Type=HW ForceField.Charge=0.417
        H 5.369672396356584 1.083579289595597 13.29812675901156 ForceField.Type=HW ForceField.Charge=0.417
        O 5.466225740442832 5.993496838594383 3.826556322251799 ForceField.Type=OW ForceField.Charge=-0.834
        H 4.698727888102344 0.3975271277797332 3.345862980059084 ForceField.Type=HW ForceField.Charge=0.417
        H 5.678590795097077 0.7684784866437412 4.409847532491263 ForceField.Type=HW ForceField.Charge=0.417
        O 5.308923379504553 4.019057754213824 12.30890772203438 ForceField.Type=OW ForceField.Charge=-0.834
        H 4.947834409182614 3.130476283899112 12.54309710913415 ForceField.Type=HW ForceField.Charge=0.417
        H 4.507522419442757 4.539019535791192 12.2645840765018 ForceField.Type=HW ForceField.Charge=0.417
        O 0.5246189538877805 2.11444774753437 5.633847854203101 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.4798849442749721 2.981151461335371 6.079301327936341 ForceField.Type=HW ForceField.Charge=0.417
        H 0.1746402018766592 1.514431334250879 6.300570842670576 ForceField.Type=HW ForceField.Charge=0.417
        O 0.7532540279197091 1.724451179192218 15.61527324613631 ForceField.Type=OW ForceField.Charge=-0.834
        H 1.102056945497837 1.316928430428979 16.4265184790613 ForceField.Type=HW ForceField.Charge=0.417
        H 0.7972987223637659 2.664620455755369 15.82837324135009 ForceField.Type=HW ForceField.Charge=0.417
        O 0.9018681569836331 3.656882816889424 3.402291999454496 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.4026827408969129 4.459281159084593 3.630496065021554 ForceField.Type=HW ForceField.Charge=0.417
        H 0.5526952071471827 3.010307298769798 4.02681327048169 ForceField.Type=HW ForceField.Charge=0.417
        O 4.863403687774274 0.9028609469360367 0.558021938995828 ForceField.Type=OW ForceField.Charge=-0.834
        H 3.996919626672331 1.319885981099747 0.7021352725814224 ForceField.Type=HW ForceField.Charge=0.417
        H 4.689645967466681 5.893899093581502 0.6975454946677974 ForceField.Type=HW ForceField.Charge=0.417
        O 0.6125422333973695 4.471658272146462 15.20199530885226 ForceField.Type=OW ForceField.Charge=-0.834
        H 5.757700386271314 4.628928071051244 14.72292292770029 ForceField.Type=HW ForceField.Charge=0.417
        H 0.4490343750054829 4.454474545873454 16.17738658790725 ForceField.Type=HW ForceField.Charge=0.417
        O 0.6366518329200573 4.752909677710202 9.952302751289087 ForceField.Type=OW ForceField.Charge=-0.834
        H 0.9300822912913463 5.41930018832904 10.60529619335514 ForceField.Type=HW ForceField.Charge=0.417
        H 0.09536209451509015 4.185550544321785 10.53226705454537 ForceField.Type=HW ForceField.Charge=0.417
        O 1.973521991893519 4.136851000321193 0.9323878175561867 ForceField.Type=OW ForceField.Charge=-0.834
        H 1.958801293505274 4.982561965771234 0.471890391343053 ForceField.Type=HW ForceField.Charge=0.417
        H 1.500116599864132 4.339560036883479 1.75645305944565 ForceField.Type=HW ForceField.Charge=0.417
        O 4.912922003545703 4.048464241294365 1.122473368577996 ForceField.Type=OW ForceField.Charge=-0.834
        H 5.258380821455368 3.682506383808921 1.973193063451997 ForceField.Type=HW ForceField.Charge=0.417
        H 3.961611234368922 3.899191033445874 1.229424459523123 ForceField.Type=HW ForceField.Charge=0.417
        O 3.68652628868048 1.595323244715648 10.09294194440881 ForceField.Type=OW ForceField.Charge=-0.834
        H 2.716738755606992 1.65284789480996 9.785054765941782 ForceField.Type=HW ForceField.Charge=0.417
        H 3.717820588318913 1.493830202223642 11.06393026159224 ForceField.Type=HW ForceField.Charge=0.417
        O 2.322791158566589 4.525494429167161 13.13712396271713 ForceField.Type=OW ForceField.Charge=-0.834
        H 1.684244052546178 4.292584907055148 12.46669482040438 ForceField.Type=HW ForceField.Charge=0.417
        H 1.713022179274242 4.656194947137037 13.91280050486365 ForceField.Type=HW ForceField.Charge=0.417
        O 1.531573866397418 1.053168936091418 0.3718291838441352 ForceField.Type=OW ForceField.Charge=-0.834
        H 1.522293261003647 2.011561261078512 0.5486474226338094 ForceField.Type=HW ForceField.Charge=0.417
        H 0.5646086101553838 0.8818604575421908 0.3306918436047358 ForceField.Type=HW ForceField.Charge=0.417
        O 3.418411885135323 2.370515151746976 3.6606774188027 ForceField.Type=OW ForceField.Charge=-0.834
        H 2.581438388703135 2.63052383259613 3.242613350177665 ForceField.Type=HW ForceField.Charge=0.417
        H 3.469669459117304 3.023705121589598 4.375706151452582 ForceField.Type=HW ForceField.Charge=0.417
    End
    Lattice
        5.952918007318143 0.0 0.0
        0.0 5.952918007318143 0.0
        0.0 0.0 17.85875402195457
    End
    BondOrders
         1 2 1.0
         1 3 1.0
         4 5 1.0
         4 6 1.0
         7 8 1.0
         7 9 1.0
         10 11 1.0
         10 12 1.0
         13 14 1.0
         13 15 1.0
         16 17 1.0
         16 18 1.0
         19 20 1.0
         19 21 1.0 -1 0 0
         22 23 1.0 0 1 0
         22 24 1.0 0 1 0
         25 26 1.0
         25 27 1.0
         28 29 1.0
         28 30 1.0
         31 32 1.0
         31 33 1.0
         34 35 1.0
         34 36 1.0
         37 38 1.0
         37 39 1.0 0 -1 0
         40 41 1.0 -1 0 0
         40 42 1.0
         43 44 1.0
         43 45 1.0
         46 47 1.0
         46 48 1.0
         49 50 1.0
         49 51 1.0
         52 53 1.0
         52 54 1.0
         55 56 1.0
         55 57 1.0
         58 59 1.0
         58 60 1.0
         61 62 1.0
         61 63 1.0
    End
End

Engine ForceField
    Type Amber95
EndEngine
eor

# The following script extracts the induced velocity amplitude and calculates the corresponding viscosity
#
# The --test argument is only useful for regression testing
"$AMSBIN/amspython" "${TEST_DIRECTORY:-.}/analyze_cosine_shear.py" --test --from 41 --block 100 ams.results/
echo "END OF TEST"