[Swift-commit] r5893 - in SwiftApps/CMTS/dimers: . bin conf input_files src
davidk at ci.uchicago.edu
davidk at ci.uchicago.edu
Wed Aug 8 10:13:06 CDT 2012
Author: davidk
Date: 2012-08-08 10:12:55 -0500 (Wed, 08 Aug 2012)
New Revision: 5893
Added:
SwiftApps/CMTS/dimers/README
SwiftApps/CMTS/dimers/TODO
SwiftApps/CMTS/dimers/bin/
SwiftApps/CMTS/dimers/bin/insert_molecules
SwiftApps/CMTS/dimers/bin/restart2data
SwiftApps/CMTS/dimers/conf/
SwiftApps/CMTS/dimers/conf/beagle.cf
SwiftApps/CMTS/dimers/conf/beagle.xml
SwiftApps/CMTS/dimers/conf/grotthuss-ssh.cf
SwiftApps/CMTS/dimers/conf/grotthuss-ssh.xml
SwiftApps/CMTS/dimers/conf/grotthuss.cf
SwiftApps/CMTS/dimers/conf/grotthuss.xml
SwiftApps/CMTS/dimers/conf/local.cf
SwiftApps/CMTS/dimers/conf/local.xml
SwiftApps/CMTS/dimers/conf/makena-ssh.cf
SwiftApps/CMTS/dimers/conf/makena-ssh.xml
SwiftApps/CMTS/dimers/conf/makena.cf
SwiftApps/CMTS/dimers/conf/makena.xml
SwiftApps/CMTS/dimers/conf/makgroth.cf
SwiftApps/CMTS/dimers/conf/makgroth.xml
SwiftApps/CMTS/dimers/conf/mcs.cf
SwiftApps/CMTS/dimers/conf/mcs.conf
SwiftApps/CMTS/dimers/conf/pads-ssh.cf
SwiftApps/CMTS/dimers/conf/pads-ssh.xml
SwiftApps/CMTS/dimers/conf/pads.cf
SwiftApps/CMTS/dimers/conf/pads.xml
SwiftApps/CMTS/dimers/dimers.swift
SwiftApps/CMTS/dimers/input_files/
SwiftApps/CMTS/dimers/input_files/end_20.restart
SwiftApps/CMTS/dimers/input_files/end_25.restart
SwiftApps/CMTS/dimers/input_files/end_30.restart
SwiftApps/CMTS/dimers/input_files/end_35.restart
SwiftApps/CMTS/dimers/input_files/end_40.restart
SwiftApps/CMTS/dimers/input_files/end_45.restart
SwiftApps/CMTS/dimers/input_files/end_50.restart
SwiftApps/CMTS/dimers/input_files/end_55.restart
SwiftApps/CMTS/dimers/input_files/end_60.restart
SwiftApps/CMTS/dimers/input_files/end_65.restart
SwiftApps/CMTS/dimers/input_files/end_70.restart
SwiftApps/CMTS/dimers/rundimers.sh
SwiftApps/CMTS/dimers/src/
SwiftApps/CMTS/dimers/src/GeometryUtil.cpp
SwiftApps/CMTS/dimers/src/GeometryUtil.h
SwiftApps/CMTS/dimers/src/Makefile
SwiftApps/CMTS/dimers/src/Ran1.cpp
SwiftApps/CMTS/dimers/src/Ran1.h
SwiftApps/CMTS/dimers/src/StringUtil.cpp
SwiftApps/CMTS/dimers/src/StringUtil.h
SwiftApps/CMTS/dimers/src/definitions.h
SwiftApps/CMTS/dimers/src/insert_molecules
SwiftApps/CMTS/dimers/src/insert_molecules.cpp
SwiftApps/CMTS/dimers/src/restart2data
SwiftApps/CMTS/dimers/src/restart2data.cpp
Log:
Initial commit of general directory layouts, source files, and config files
Added: SwiftApps/CMTS/dimers/README
===================================================================
Added: SwiftApps/CMTS/dimers/TODO
===================================================================
Added: SwiftApps/CMTS/dimers/bin/insert_molecules
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/bin/insert_molecules
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/bin/restart2data
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/bin/restart2data
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/conf/beagle.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/beagle.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/beagle.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,10 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=3
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+
+#app sumloss=$PWD/../sumloss.sh
+#app evolve=$PWD/../evolve.sh
Added: SwiftApps/CMTS/dimers/conf/beagle.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/beagle.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/beagle.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,23 @@
+<config>
+ <pool handle="beagle">
+ <execution provider="coaster" jobmanager="local:pbs"/>
+
+ <profile namespace="env" key="SWIFT_GEN_SCRIPTS">KEEP</profile>
+
+ <profile namespace="globus" key="jobsPerNode">1</profile>
+ <profile namespace="globus" key="lowOverAllocation">100</profile>
+ <profile namespace="globus" key="highOverAllocation">100</profile>
+ <profile namespace="globus" key="providerAttributes">pbs.aprun;pbs.mpp;depth=24</profile>
+ <profile namespace="globus" key="maxTime">10000</profile>
+ <profile namespace="globus" key="maxWallTime">01:30:00</profile>
+ <profile namespace="globus" key="slots">50</profile>
+ <profile namespace="globus" key="nodeGranularity">2</profile>
+ <profile namespace="globus" key="maxNodes">2</profile>
+ <profile namespace="globus" key="queue">route</profile>
+ <profile namespace="karajan" key="jobThrottle">9.59</profile>
+ <profile namespace="karajan" key="initialScore">10000</profile>
+ <filesystem provider="local"/>
+ <workdirectory>_WORK_/beagle</workdirectory>
+ </pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/conf/grotthuss-ssh.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/grotthuss-ssh.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/grotthuss-ssh.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=false
+use.provider.staging=true
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+tcp.port.range=5000,51000
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/grotthuss-ssh.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/grotthuss-ssh.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/grotthuss-ssh.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,20 @@
+<config>
+<pool handle="grotthuss-ssh">
+ <execution provider="coaster" url="grotthuss.uchicago.edu" jobmanager="ssh-cl:pbs"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="HighOverAllocation">100</profile>
+ <profile namespace="globus" key="LowOverAllocation">100</profile>
+ <profile namespace="globus" key="maxWallTime">00:50:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">1.01</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">1</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
Added: SwiftApps/CMTS/dimers/conf/grotthuss.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/grotthuss.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/grotthuss.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/grotthuss.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/grotthuss.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/grotthuss.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,19 @@
+<config>
+<pool handle="grotthuss">
+ <execution jobmanager="local:pbs" provider="pbs" url="none"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="maxWallTime">01:00:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">8</profile>
+ <profile key="slots" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">1</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">5.99</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">1</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/conf/local.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/local.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/local.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,10 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+
+#app rmsd=$PWD/../rmsd.sh
Added: SwiftApps/CMTS/dimers/conf/local.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/local.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/local.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,14 @@
+<config>
+ <pool handle="localhost">
+ <filesystem provider="local" />
+ <execution provider="coaster" jobmanager="local:local"/>
+ <profile namespace="karajan" key="jobthrottle">2.55</profile>
+ <profile namespace="karajan" key="initialScore">10000</profile>
+ <profile namespace="globus" key="jobsPerNode">4</profile>
+ <profile namespace="globus" key="slots">8</profile>
+ <profile namespace="globus" key="maxTime">1000</profile>
+ <profile namespace="globus" key="nodeGranularity">1</profile>
+ <profile namespace="globus" key="maxNodes">4</profile>
+ <workdirectory>/tmp</workdirectory>
+ </pool>
+</config>
Added: SwiftApps/CMTS/dimers/conf/makena-ssh.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/makena-ssh.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makena-ssh.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=false
+use.provider.staging=true
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+tcp.port.range=5000,51000
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/makena-ssh.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/makena-ssh.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makena-ssh.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,21 @@
+<config>
+<pool handle="makena-ssh">
+ <execution provider="coaster" url="makena.uchicago.edu" jobmanager="ssh:pbs"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="HighOverAllocation">100</profile>
+ <profile namespace="globus" key="LowOverAllocation">100</profile>
+ <profile namespace="globus" key="maxWallTime">00:50:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">1.01</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">2</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/conf/makena.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/makena.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makena.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/makena.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/makena.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makena.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,20 @@
+<config>
+<pool handle="makena">
+ <execution jobmanager="local:pbs" provider="pbs" url="none"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="maxWallTime">01:00:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">8</profile>
+ <profile key="slots" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">2</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">5.99</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">2</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/conf/makgroth.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/makgroth.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makgroth.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=0
+lazy.errors=false
+use.provider.staging=true
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+tcp.port.range=5000,51000
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/makgroth.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/makgroth.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/makgroth.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,41 @@
+<config>
+
+<pool handle="makena-ssh">
+ <execution provider="coaster" url="makena.uchicago.edu" jobmanager="ssh-cl:pbs"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="HighOverAllocation">100</profile>
+ <profile namespace="globus" key="LowOverAllocation">100</profile>
+ <profile namespace="globus" key="maxWallTime">00:50:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">1.01</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">2</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+
+<pool handle="grotthuss-ssh">
+ <execution provider="coaster" url="grotthuss.uchicago.edu" jobmanager="ssh-cl:pbs"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="HighOverAllocation">100</profile>
+ <profile namespace="globus" key="LowOverAllocation">100</profile>
+ <profile namespace="globus" key="maxWallTime">00:50:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">devel</profile>
+ <profile key="jobThrottle" namespace="karajan">1.01</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">1</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+
+</config>
Added: SwiftApps/CMTS/dimers/conf/mcs.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/mcs.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/mcs.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,7 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=3
+lazy.errors=true
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
Added: SwiftApps/CMTS/dimers/conf/mcs.conf
===================================================================
--- SwiftApps/CMTS/dimers/conf/mcs.conf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/mcs.conf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,45 @@
+# Keep all interesting settings in one place
+# User should modify this to fit environment
+
+# Location of SWIFT. If empty, PATH is referenced
+export SWIFT=
+
+# Where to place/launch worker.pl on the remote machine for sites.xml
+export WORKER_LOCATION=_RUNDIR_
+
+# How to launch workers: local, ssh, or cobalt
+export WORKER_MODE=ssh
+
+# Worker logging setting passed to worker.pl for sites.xml
+export WORKER_LOGGING_LEVEL=INFO
+
+# User name to use for all systems
+export WORKER_USERNAME=$USER
+
+# Worker host names for ssh
+export WORKER_HOSTS="crush thwomp stomp crank steamroller grind churn trounce thrash vanquish"
+
+# Directory to keep log files, relative to working directory when launching start-coaster-service
+export LOG_DIR=logs
+export WORKER_LOG_DIR=_RUNDIR_
+
+# Manually define ports. If not specified, ports will be automatically generated
+export LOCAL_PORT=
+export SERVICE_PORT=
+
+# Set shared filesystem to no since work will be done in local /sandbox directory
+export SHARED_FILESYSTEM=yes
+
+# start-coaster-service tries to automatically detect IP address.
+# Specify here if auto detection is not working correctly
+export IPADDR=
+
+# Below are various settings to give information about how to create sites.xml
+#export WORK=_RUNDIR_
+export JOBS_PER_NODE=1
+# export JOB_THROTTLE=$( echo "scale=5; $( echo $WORKER_HOSTS | wc -w )/100 - 0.00001"|bc )
+
+# Swift applications
+#app sumloss=$PWD/../sumloss.sh
+#app evolve=$PWD/../evolve.sh
+
Added: SwiftApps/CMTS/dimers/conf/pads-ssh.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/pads-ssh.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/pads-ssh.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=10
+lazy.errors=false
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/pads-ssh.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/pads-ssh.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/pads-ssh.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,21 @@
+<config>
+<pool handle="makena-ssh">
+ <execution provider="coaster" url="login1.pads.ci.uchicago.edu" jobmanager="ssh-cl:pbs"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="HighOverAllocation">100</profile>
+ <profile namespace="globus" key="LowOverAllocation">100</profile>
+ <profile namespace="globus" key="maxWallTime">00:50:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">1</profile>
+ <profile key="maxNodes" namespace="globus">2</profile>
+ <profile key="queue" namespace="globus">fast</profile>
+ <profile key="jobThrottle" namespace="karajan">1.01</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">2</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/conf/pads.cf
===================================================================
--- SwiftApps/CMTS/dimers/conf/pads.cf (rev 0)
+++ SwiftApps/CMTS/dimers/conf/pads.cf 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,11 @@
+wrapperlog.always.transfer=true
+sitedir.keep=true
+execution.retries=10
+lazy.errors=false
+status.mode=provider
+use.provider.staging=false
+provider.staging.pin.swiftfiles=false
+use.wrapper.staging=false
+
+#app rmsd=$PWD/../rmsd.sh
+#app plot_pd=$PWD/../plot_pd.pl
Added: SwiftApps/CMTS/dimers/conf/pads.xml
===================================================================
--- SwiftApps/CMTS/dimers/conf/pads.xml (rev 0)
+++ SwiftApps/CMTS/dimers/conf/pads.xml 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,20 @@
+<config>
+<pool handle="makena">
+ <execution jobmanager="local:pbs" provider="pbs" url="none"/>
+ <filesystem provider="local" url="none" />
+ <profile namespace="globus" key="maxWallTime">01:00:00</profile>
+ <profile namespace="globus" key="maxTime">3600</profile>
+ <profile key="jobsPerNode" namespace="globus">8</profile>
+ <profile key="slots" namespace="globus">1</profile>
+ <profile key="nodeGranularity" namespace="globus">2</profile>
+ <profile key="maxNodes" namespace="globus">1</profile>
+ <profile key="queue" namespace="globus">fast</profile>
+ <profile key="jobThrottle" namespace="karajan">5.99</profile>
+ <profile key="initialScore" namespace="karajan">10000</profile>
+ <profile key="count" namespace="globus">1</profile>
+ <profile key="jobType" namespace="globus">single</profile>
+ <profile key="ppn" namespace="globus">8</profile>
+ <workdirectory>_WORK_</workdirectory>
+</pool>
+</config>
+
Added: SwiftApps/CMTS/dimers/dimers.swift
===================================================================
--- SwiftApps/CMTS/dimers/dimers.swift (rev 0)
+++ SwiftApps/CMTS/dimers/dimers.swift 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,13 @@
+type file;
+
+app insert_molecules
+
+app restart2data
+
+app lmp_john
+
+int start_ndimers=20;
+int max_ndimers=75;
+int delta_ndimers=5;
+int virion_radius=300;
+
Added: SwiftApps/CMTS/dimers/input_files/end_20.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_20.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_25.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_25.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_30.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_30.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_35.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_35.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_40.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_40.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_45.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_45.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_50.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_50.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_55.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_55.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_60.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_60.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_65.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_65.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/input_files/end_70.restart
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/input_files/end_70.restart
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/rundimers.sh
===================================================================
--- SwiftApps/CMTS/dimers/rundimers.sh (rev 0)
+++ SwiftApps/CMTS/dimers/rundimers.sh 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,90 @@
+#!/bin/bash
+
+# Usage: swiftopt.sh [-s sitename] [-w]
+#
+
+usage="$0 [-s sitename] [-p paramfile]"
+
+# Function to run Swift
+runswift() {
+ swift -tc.file tc.data -sites.file $1 -config cf dimers.swift 2>&1 | tee swift.out
+}
+
+# Default settings
+execsite=local
+paramfile=local
+ram=2048M
+
+# Process command line arguments
+while [ $# -gt 0 ]; do
+ case $1 in
+ -s) execsite=$2; shift 2;;
+ *) echo $usage 1>&2
+ exit 1;;
+ esac
+done
+
+# Create next unique run id and run directory
+rundir=$( echo run??? | sed -e 's/^.*run//' | awk '{ printf("run%03d\n", $1+1)}' )
+
+# Exit if rundir already exits. Something is funky
+if [ -d $rundir ];
+then
+ echo "$rundir already exists! exiting." >&2
+ exit 2
+else
+ mkdir $rundir
+fi
+
+# Copy input files
+cp input_files/* $rundir
+echo Run directory $rundir: site=$execsite paramfile=$paramfile
+
+# Report an error if configuration files are missing
+if [ ! -f "conf/$execsite.xml" ] && [ ! -f "conf/$execsite.conf" ]; then
+ echo Unable to find requested configuration file for site $execsite
+ exit 1
+fi
+
+# Use start-coaster-service if site is a .conf file
+if [ -f "conf/$execsite.conf" ]; then
+ USE_SCS=1
+fi
+
+# Check for missing .cf files
+if [ -f "conf/$execsite.xml" ] && [ ! -f "conf/$execsite.cf" ]; then
+ echo Missing configuration file $execsite.cf
+fi
+
+# Do the run
+cd $rundir
+cp ../rmsd.swift .
+export WORK=$PWD/swiftwork
+mkdir -p $PWD/swiftwork/workers
+
+# Use start-coaster-service if the site uses a .conf file
+if [ "$USE_SCS" == "1" ]; then
+ cp ../conf/$execsite.conf coaster-service.conf
+ cp ../conf/$execsite.cf cf
+ sed -i -e "s at _RUNDIR_@$rundir@" coaster-service.conf
+ start-coaster-service
+fi
+
+# Run gensites
+if [ ! "$USE_SCS" == 1 ]; then
+ cp ../conf/$execsite.cf cf
+ gensites -p ../conf/$execsite.cf ../conf/$execsite.xml > $execsite.xml
+fi
+
+echo "Run dir=$rundir" >> ABOUT
+echo "Work dir=$WORK" >> ABOUT
+echo "Total jobs=$total_jobs" >> ABOUT
+
+if [ "$USE_SCS" == "1" ]; then
+ runswift "sites.xml"
+ stop-coaster-service
+else
+ runswift "$execsite.xml"
+fi
+
+exit
Property changes on: SwiftApps/CMTS/dimers/rundimers.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/GeometryUtil.cpp
===================================================================
--- SwiftApps/CMTS/dimers/src/GeometryUtil.cpp (rev 0)
+++ SwiftApps/CMTS/dimers/src/GeometryUtil.cpp 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,721 @@
+#include "GeometryUtil.h"
+
+/*
+ Here theta is the "horizontal" resolution, and phi is the "vertical resolution"
+ We add single point "caps" on the north and south poles of the sphere - these are included
+ in the phi resolution, so the number of vertical circles which for the sphere is phi_res-2.
+*/
+FLOAT_TYPE * GenerateSphere( int theta_res, int phi_res, FLOAT_TYPE sphere_radius, int *npoints )
+{
+ double theta, phi, delta_phi;
+ double circle_radius, pi, twopi, x, y, z;
+ FLOAT_TYPE *points; // FLOAT_TYPE, as we pass it back.
+ int i, j, k;
+
+ assert( phi_res > 2 ); // ensure we have at least one genuine circle on the vertical axis!
+
+ pi = M_PI;
+ twopi = 2.0*M_PI;
+
+ phi_res -= 2;
+
+ *npoints = 2 + (phi_res*theta_res);
+ points = (FLOAT_TYPE *) malloc( sizeof(FLOAT_TYPE)*(*npoints)*3 );
+
+ k = 0;
+
+ // add north pole - single point
+ points[k*3 +0] = 0.0;
+ points[k*3 +1] = 0.0;
+ points[k*3 +2] = sphere_radius;
+ k++;
+
+ // add list of circles on vertical axis
+ delta_phi = pi/(phi_res+1);
+ for( i=0; i<phi_res; i++ )
+ {
+ phi = delta_phi*(i+1);
+ z = sphere_radius*cos( phi );
+ circle_radius = sphere_radius*sin( phi );
+
+ for( j=0; j<theta_res; j++ )
+ {
+ theta = (twopi/theta_res)*j;
+ x = circle_radius*sin(theta);
+ y = circle_radius*cos(theta);
+
+ points[k*3 +0] = x;
+ points[k*3 +1] = y;
+ points[k*3 +2] = z;
+ k++;
+ }
+ }
+
+ // add south pole - single point
+ points[k*3 +0] = 0.0;
+ points[k*3 +1] = 0.0;
+ points[k*3 +2] = -sphere_radius;
+ k++;
+
+ assert( k == *npoints );
+
+ return points;
+}
+
+
+void CentreOfMass( int N, FLOAT_TYPE *xyz, FLOAT_TYPE *mass, FLOAT_TYPE *COM )
+{
+ int i;
+ double m, total_mass;
+
+ assert( N > 0 );
+ assert( xyz != NULL );
+ assert( COM != NULL );
+
+ COM[0] = 0.0;
+ COM[1] = 0.0;
+ COM[2] = 0.0;
+ total_mass = 0.0;
+
+ if( mass == NULL )
+ {
+ // assume mass of 1.0
+ m = 1.0;
+ for( i=0; i<N; i++ )
+ {
+ COM[0] += m * xyz[i*3 +0];
+ COM[1] += m * xyz[i*3 +1];
+ COM[2] += m * xyz[i*3 +2];
+ total_mass += m;
+ }
+ }
+ else
+ {
+ for( i=0; i<N; i++ )
+ {
+ m = mass[i];
+ COM[0] += m * xyz[i*3 +0];
+ COM[1] += m * xyz[i*3 +1];
+ COM[2] += m * xyz[i*3 +2];
+ total_mass += m;
+ }
+ }
+
+ COM[0] /= total_mass;
+ COM[1] /= total_mass;
+ COM[2] /= total_mass;
+}
+
+/*
+ I want to maintain the full range of coords for theta and phi; is we simply use the common vector angle difference
+ method, we only get values of 0->pi, ie the minimum angle. For the case of theta, this gives us a circle around the zenith axis
+ but to properly define the phi we need the full range of 0->2pi so we don't have uncertainty as to which "side" of the circle
+ phi actually lies on.
+
+ This is an issue, as which "side" of the circle phi is on dependes on whether we're looking from above or below the origin of
+ the coordinate system, as that reflects phi! Therefore, we make an assumption here; specifically, that the zenith axis defines
+ a normal to the plane formed by the remaining two cartesian axes, and we use this as "up". If we then project the 3D vector
+ onto a 2D plane defined by the two potential azimuthal axes (by killing the zenith axis component), we can actually get the
+ full range of 0->2pi given that "up" is the zenith axis, and IT IS SIGNED!
+
+ Here we use zenith axis = z and azimuthal axis = x for convenience and easy interpretation.
+
+ "Conventional" method (phi only on [0,pi]):
+
+ x,y,z = v2-v1
+ r = sqrt( dot(v1,v2) )
+ theta = acos( z / r )
+ phi = atan( y, x )
+*/
+void CartesianToSpherical( FLOAT_TYPE *p1, FLOAT_TYPE *p2, FLOAT_TYPE *r, FLOAT_TYPE *theta, FLOAT_TYPE *phi )
+{
+ double vec[3], axis[3], r2, x, y;
+ double min_r = 1e-20;
+
+ assert( p1 != NULL );
+ assert( p2 != NULL );
+ assert( r != NULL );
+ assert( theta != NULL );
+ assert( phi != NULL );
+
+ vec[0] = p2[0] - p1[0];
+ vec[1] = p2[1] - p1[1];
+ vec[2] = p2[2] - p1[2];
+
+ *r = sqrt( DOT_PRODUCT(vec,vec) );
+ // is the 3D point too close to the origin?
+ if( *r < min_r )
+ {
+ *theta = 0.0;
+ *phi = 0.0;
+ }
+ else
+ {
+ vec[0] /= *r;
+ vec[1] /= *r;
+ vec[2] /= *r;
+
+ // theta is the angle between the vector and the zenith axis - this is easy.
+ axis[0] = axis[1] = axis[2] = 0.0;
+ axis[2] = 1.0; // zenith axis is z
+ *theta = acos( DOT_PRODUCT(axis,vec) );
+
+ // now we get clever; we project the vector onto the 2D plane by killing the zenith component. Obviously,
+ // this vector may no longer be a unit vector, so we renormalise. We can then use the atan2 function to get
+ // the angle in the range -pi->pi (atan2 is also more accurate than other methods)
+ axis[0] = axis[1] = axis[2] = 0.0;
+ axis[0] = 1.0; // azimuthal axis is x
+
+ vec[2] = 0.0; // kill the zenith component of the vector to project onto azimuthal plane
+ r2 = sqrt( DOT_PRODUCT(vec,vec) ); // renomalise value
+ // is the projected point too close to the 2D origin?
+ if( r2 < min_r ) *phi = 0.0;
+ else
+ {
+ // use the atan2 function, with x and y components!
+ vec[0] /= r2;
+ vec[1] /= r2;
+ vec[2] /= r2;
+ //
+ x = vec[0];
+ y = vec[1];
+ *phi = atan2( x, y );
+ //*phi = acos( axis[0]*vec[0] + axis[1]*vec[1] + axis[2]*vec[2] ); // old version; only 0 -> pi, so merges data symmetrically.
+ }
+ }
+}
+void SphericalToCartesian( FLOAT_TYPE r, FLOAT_TYPE theta, FLOAT_TYPE phi, FLOAT_TYPE *x, FLOAT_TYPE *y, FLOAT_TYPE *z )
+{
+ assert( x != NULL );
+ assert( y != NULL );
+ assert( z != NULL );
+
+ // deffo right! unaffected by phi calculations.
+ *z = r * cos(theta);
+
+ // but it fucking beats me why the x and y components have switched here vs "standard" spherical. :/
+ *x = r * sin(theta)*sin( phi );
+ *y = r * sin(theta)*cos( phi );
+
+ // the above actually avoids singularities that catch out the standard spherical coords method,
+ // so we don't need to check for any weirdness as we have phi defined with no weirdness.
+
+ // alternatively ...
+ //FLOAT_TYPE u;
+ //u = sqrt( r*r - (*z)*(*z) );
+ //*x = u * sin( phi );
+ //*y = u * cos( phi );
+}
+
+/*
+// conventional Cartesian => spherical coords conversion
+void _CartesianToSpherical( FLOAT_TYPE *p1, FLOAT_TYPE *p2, FLOAT_TYPE *r, FLOAT_TYPE *theta, FLOAT_TYPE *phi )
+{
+ double x, y, z;
+
+ x = p2[0] - p1[0];
+ y = p2[1] - p1[1];
+ z = p2[2] - p1[2];
+
+ *r = sqrt( x*x + y*y + z*z );
+ *theta = acos( z / *r );
+ *phi = atan( y/x );
+}
+// conventional spherical coords => Cartesian conversion
+void _SphericalToCartesian( FLOAT_TYPE r, FLOAT_TYPE theta, FLOAT_TYPE phi, FLOAT_TYPE *x, FLOAT_TYPE *y, FLOAT_TYPE *z )
+{
+ *x = r * sin(theta)*cos(phi);
+ *y = r * sin(theta)*sin(phi);
+ *z = r * cos(theta);
+}
+// check we recover correct coordinates vs original, and against standard spherical methods marked with underscore.
+// note the problems in the original method that we avoid around the origin!
+void test_spherical( int dim )
+{
+ int x, y, z;
+ FLOAT_TYPE s_r, s_t, s_p; // need to be FLOAT_TYPE, as we're passing into and out of the FLOAT_TYPE routines.
+ FLOAT_TYPE v1[3], v2[3]; // as above.
+
+ v1[0] = v1[1] = v1[2] = 0.0;
+
+ for( x=-dim; x<=dim; x++ )
+ for( y=-dim; y<=dim; y++ )
+ for( z=-dim; z<=dim; z++ )
+ {
+ v2[0] = x;
+ v2[1] = y;
+ v2[2] = z;
+ printf( "%+.12f %+.12f %+.12f\n", v2[0], v2[1], v2[2] );
+
+ _CartesianToSpherical( v1, v2, &s_r, &s_t, &s_p );
+ printf( "\t %+.6f %+.6f %+.6f => ", s_r, s_t, s_p );
+ _SphericalToCartesian( s_r, s_t, s_p, &v2[0], &v2[1], &v2[2] );
+ printf( "%+.12f %+.12f %+.12f\n", v2[0], v2[1], v2[2] );
+
+ v2[0] = x;
+ v2[1] = y;
+ v2[2] = z;
+ CartesianToSpherical( v1, v2, &s_r, &s_t, &s_p );
+ printf( "\t %+.6f %+.6f %+.6f => ", s_r, s_t, s_p );
+ SphericalToCartesian( s_r, s_t, s_p, &v2[0], &v2[1], &v2[2] );
+ printf( "%+.12f %+.12f %+.12f\n", v2[0], v2[1], v2[2] );
+ }
+}
+*/
+
+void Translate( int N, FLOAT_TYPE *xyz, FLOAT_TYPE *offset )
+{
+ int i;
+
+ assert( N > 0 );
+ assert( xyz != NULL );
+ assert( offset != NULL );
+
+ for( i=0; i<N; i++ )
+ {
+ xyz[i*3 +0] += offset[0];
+ xyz[i*3 +1] += offset[1];
+ xyz[i*3 +2] += offset[2];
+ }
+}
+// pass point = NULL for rotation about system origin, else rotation axis will pass through that point.
+void RotateAboutPoint( int N, FLOAT_TYPE *xyz, FLOAT_TYPE *point, FLOAT_TYPE *axis, FLOAT_TYPE theta )
+{
+ int i;
+ double a, b, c, x, y, z, rx, ry, rz, u, v, w, u2, v2, w2, K, L, ct, st;
+
+ assert( N > 0 );
+ assert( xyz != NULL );
+ assert( axis != NULL );
+
+ a = 0.0;
+ b = 0.0;
+ c = 0.0;
+
+ if( point != NULL )
+ {
+ a = point[0];
+ b = point[1];
+ c = point[2];
+ }
+
+ ct = cos( theta );
+ st = sin( theta );
+
+ u = axis[0];
+ v = axis[1];
+ w = axis[2];
+ u2 = u*u;
+ v2 = v*v;
+ w2 = w*w;
+ K = u2 + v2 + w2;
+ L = sqrt( K );
+
+ for( i=0; i<N; i++ )
+ {
+ x = xyz[i*3 +0] - a;
+ y = xyz[i*3 +1] - b;
+ z = xyz[i*3 +2] - c;
+
+ rx = a*(v2+w2) + u*(-b*v -c*w + u*x + v*y + w*z) + ( (x-a)*(v2+w2) + u*(b*v + c*w - v*y - w*z) ) * ct + L*(b*w - c*v - w*y + v*z)*st;
+ rx = rx / K;
+
+ ry = b*(u2+w2) + v*(-a*u -c*w + u*x + v*y + w*z) + ( (y-b)*(u2+w2) + v*(a*u + c*w - u*x - w*z) ) * ct + L*(-a*w + c*u + w*x - u*z)*st;
+ ry = ry / K;
+
+ rz = c*(u2+v2) + w*(-a*u -b*v + u*x + v*y + w*z) + ( (z-c)*(u2+v2) + w*(a*u + b*v - u*x - v*y) ) * ct + L*(a*v - b*u - v*x + u*y)*st;
+ rz = rz / K;
+
+ xyz[i*3 +0] = rx + a;
+ xyz[i*3 +1] = ry + b;
+ xyz[i*3 +2] = rz + c;
+ }
+}
+
+
+
+FLOAT_TYPE l2_norm_difference( int len, FLOAT_TYPE *v1, FLOAT_TYPE *v2 )
+{
+ int i;
+ double d, acc;
+
+ acc = 0.0;
+ for( i=0; i<len; i++ )
+ {
+ d = v1[i] - v2[i];
+ acc += d*d;
+ }
+ return sqrt( acc );
+}
+
+
+void Subdivider::InitTetrahedron()
+{
+ float sqrt3 = 1 / sqrt(3.0);
+ float tetrahedron_vertices[] = {
+ sqrt3, sqrt3, sqrt3,
+ -sqrt3, -sqrt3, sqrt3,
+ -sqrt3, sqrt3, -sqrt3,
+ sqrt3, -sqrt3, -sqrt3 };
+ int tetrahedron_faces[] = {0, 2, 1, 0, 1, 3, 2, 3, 1, 3, 2, 0};
+
+ n_vertices = 4;
+ n_faces = 4;
+ n_edges = 6;
+ vertices = (float*) malloc( 3*n_vertices*sizeof(float) );
+ faces = (int*) malloc( 3*n_faces*sizeof(int) );
+ memcpy( (void*) vertices, (void*) tetrahedron_vertices, 3*n_vertices*sizeof(float) );
+ memcpy( (void*) faces, (void*) tetrahedron_faces, 3*n_faces*sizeof(int) );
+}
+void Subdivider::InitOctahedron()
+{
+ float octahedron_vertices[] = {
+ 0.0, 0.0, -1.0,
+ 1.0, 0.0, 0.0,
+ 0.0, -1.0, 0.0,
+ -1.0, 0.0, 0.0,
+ 0.0, 1.0, 0.0,
+ 0.0, 0.0, 1.0 };
+ int octahedron_faces[] = {0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1, 5, 2, 1, 5, 3, 2, 5, 4, 3, 5, 1, 4};
+
+ n_vertices = 6;
+ n_faces = 8;
+ n_edges = 12;
+ vertices = (float*) malloc( 3*n_vertices*sizeof(float) );
+ faces = (int*) malloc( 3*n_faces*sizeof(int) );
+ memcpy( (void*) vertices, (void*) octahedron_vertices, 3*n_vertices*sizeof(float) );
+ memcpy( (void*) faces, (void*) octahedron_faces, 3*n_faces*sizeof(int) );
+}
+void Subdivider::InitIcosahedron()
+{
+ float t = (1+sqrt(5))/2;
+ float tau = t/sqrt(1+t*t);
+ float one = 1/sqrt(1+t*t);
+
+ float icosahedron_vertices[] = {
+ tau, one, 0.0,
+ -tau, one, 0.0,
+ -tau, -one, 0.0,
+ tau, -one, 0.0,
+ one, 0.0 , tau,
+ one, 0.0 , -tau,
+ -one, 0.0 , -tau,
+ -one, 0.0 , tau,
+ 0.0 , tau, one,
+ 0.0 , -tau, one,
+ 0.0 , -tau, -one,
+ 0.0 , tau, -one };
+
+ int icosahedron_faces[] = {
+ 4, 8, 7,
+ 4, 7, 9,
+ 5, 6, 11,
+ 5, 10, 6,
+ 0, 4, 3,
+ 0, 3, 5,
+ 2, 7, 1,
+ 2, 1, 6,
+ 8, 0, 11,
+ 8, 11, 1,
+ 9, 10, 3,
+ 9, 2, 10,
+ 8, 4, 0,
+ 11, 0, 5,
+ 4, 9, 3,
+ 5, 3, 10,
+ 7, 8, 1,
+ 6, 1, 11,
+ 7, 2, 9,
+ 6, 10, 2 };
+
+ n_vertices = 12;
+ n_faces = 20;
+ n_edges = 30;
+ vertices = (float*) malloc( 3*n_vertices*sizeof(float) );
+ faces = (int*) malloc( 3*n_faces*sizeof(int) );
+ memcpy( (void*) vertices, (void*) icosahedron_vertices, 3*n_vertices*sizeof(float) );
+ memcpy( (void*) faces, (void*) icosahedron_faces, 3*n_faces*sizeof(int) );
+}
+int Subdivider::SearchMidpoint( int index_start, int index_end )
+{
+ int i, res;
+ float length;
+
+ for (i=0; i<edge_walk; i++)
+ {
+ if( (start[i] == index_start && end[i] == index_end) || (start[i] == index_end && end[i] == index_start) )
+ {
+ res = midpoint[i];
+
+ /* update the arrays */
+ start[i] = start[edge_walk-1];
+ end[i] = end[edge_walk-1];
+ midpoint[i] = midpoint[edge_walk-1];
+ edge_walk--;
+
+ return res;
+ }
+ }
+
+ /* vertex not in the list, so we add it */
+ start[edge_walk] = index_start;
+ end[edge_walk] = index_end;
+ midpoint[edge_walk] = n_vertices;
+
+ /* create new vertex */
+ vertices[3*n_vertices+0] = (vertices[3*index_start+0] + vertices[3*index_end+0]) / 2.0;
+ vertices[3*n_vertices+1] = (vertices[3*index_start+1] + vertices[3*index_end+1]) / 2.0;
+ vertices[3*n_vertices+2] = (vertices[3*index_start+2] + vertices[3*index_end+2]) / 2.0;
+
+ /* normalize the new vertex */
+ length = sqrt(
+ vertices[3*n_vertices+0] * vertices[3*n_vertices+0] +
+ vertices[3*n_vertices+1] * vertices[3*n_vertices+1] +
+ vertices[3*n_vertices+2] * vertices[3*n_vertices+2] );
+ length = 1.0/length;
+ vertices[3*n_vertices+0] *= length;
+ vertices[3*n_vertices+1] *= length;
+ vertices[3*n_vertices+2] *= length;
+
+ n_vertices++;
+ edge_walk++;
+ return midpoint[edge_walk-1];
+}
+void Subdivider::Subdivide()
+{
+ int n_vertices_new = n_vertices+2*n_edges;
+ int n_faces_new = 4*n_faces;
+ int i;
+
+ int *faces_old;
+ int a, b, c, ab_midpoint, bc_midpoint, ca_midpoint;
+
+ edge_walk = 0;
+ n_edges = 2*n_vertices + 3*n_faces;
+ start = (int*) malloc( n_edges*sizeof(int) );
+ end = (int*) malloc( n_edges*sizeof(int) );
+ midpoint = (int*) malloc( n_edges*sizeof(int) );
+
+ faces_old = (int*) malloc( 3*n_faces*sizeof(int) );
+ faces_old = (int*) memcpy( (void*)faces_old, (void*)faces, 3*n_faces*sizeof(int) );
+ vertices = (float*) realloc( (void*)vertices, 3*n_vertices_new*sizeof(float) );
+ faces = (int*) realloc( (void*)faces, 3*n_faces_new*sizeof(int) );
+ n_faces_new = 0;
+
+ for (i=0; i<n_faces; i++)
+ {
+ a = faces_old[3*i+0];
+ b = faces_old[3*i+1];
+ c = faces_old[3*i+2];
+
+ ab_midpoint = SearchMidpoint( b, a );
+ bc_midpoint = SearchMidpoint( c, b );
+ ca_midpoint = SearchMidpoint( a, c );
+
+ faces[3*n_faces_new+0] = a;
+ faces[3*n_faces_new+1] = ab_midpoint;
+ faces[3*n_faces_new+2] = ca_midpoint;
+ n_faces_new++;
+ faces[3*n_faces_new+0] = ca_midpoint;
+ faces[3*n_faces_new+1] = ab_midpoint;
+ faces[3*n_faces_new+2] = bc_midpoint;
+ n_faces_new++;
+ faces[3*n_faces_new+0] = ca_midpoint;
+ faces[3*n_faces_new+1] = bc_midpoint;
+ faces[3*n_faces_new+2] = c;
+ n_faces_new++;
+ faces[3*n_faces_new+0] = ab_midpoint;
+ faces[3*n_faces_new+1] = b;
+ faces[3*n_faces_new+2] = bc_midpoint;
+ n_faces_new++;
+ }
+ n_faces = n_faces_new;
+ free( start );
+ free( end );
+ free( midpoint );
+ free( faces_old );
+}
+Subdivider::Subdivider( int init_as, int n_subdivisions )
+{
+ int i;
+
+ n_vertices = 0;
+ n_faces = 0;
+ n_edges = 0;
+ vertices = NULL;
+ faces = NULL;
+
+ edge_walk = 0;
+ start = NULL;
+ end = NULL;
+ midpoint = NULL;
+
+ switch( init_as )
+ {
+ case INIT_TETRAHEDRON:
+ InitTetrahedron();
+ break;
+
+ case INIT_OCTAHEDRON:
+ InitOctahedron();
+ break;
+
+ case INIT_ICOSAHEDRON:
+ InitIcosahedron();
+ break;
+
+ default:
+ printf( "Subdivider::%s(): unknown init type.\n", __func__ );
+ exit( -1 );
+ break;
+ }
+
+ double rv[3], r;
+ rv[0] = 0.0;
+ rv[1] = 0.0;
+ rv[2] = 0.0;
+ for( i=0; i<n_vertices; i++ )
+ {
+ rv[0] += vertices[i*3 +0];
+ rv[1] += vertices[i*3 +1];
+ rv[2] += vertices[i*3 +2];
+ }
+ rv[0] /= n_vertices;
+ rv[1] /= n_vertices;
+ rv[2] /= n_vertices;
+
+ //printf( "Subdivider::%s(): starting vertex count is %d\n", __func__, n_vertices );
+ for( i=0; i<n_subdivisions; i++ )
+ {
+ Subdivide();
+ //printf( "Subdivider::%s(): subdivision %d resulted in %d vertices\n", __func__, i+1, n_vertices );
+ }
+
+ // renormalise - redundant, but never mind.
+ for( i=0; i<n_vertices; i++ )
+ {
+ rv[0] = vertices[i*3 +0];
+ rv[1] = vertices[i*3 +1];
+ rv[2] = vertices[i*3 +2];
+ r = sqrt( rv[0]*rv[0] + rv[1]*rv[1] + rv[2]*rv[2] );
+ vertices[i*3+0] /= r;
+ vertices[i*3+1] /= r;
+ vertices[i*3+2] /= r;
+ }
+}
+Subdivider::~Subdivider()
+{
+ if( vertices != NULL ) free( vertices );
+ if( faces != NULL ) free( faces );
+}
+FLOAT_TYPE * Subdivider::BuildShell( int n_atoms, FLOAT_TYPE *atom_pos, FLOAT_TYPE r, int *n_shell_points )
+{
+ FLOAT_TYPE * shell_points;
+ double pos[3], rv[3], r2;
+ int i, j, k, good_point;
+
+ std::vector<FLOAT_TYPE> x, y, z;
+
+ // build a list of "good" vertices in the shell
+ for( i=0; i<n_atoms; i++ )
+ {
+ for( j=0; j<n_vertices; j++ )
+ {
+ pos[0] = atom_pos[i*3 +0] + vertices[j*3 +0]*r;
+ pos[1] = atom_pos[i*3 +1] + vertices[j*3 +1]*r;
+ pos[2] = atom_pos[i*3 +2] + vertices[j*3 +2]*r;
+
+ // check whether vertex collides with another atom
+ good_point = 1;
+ for( k=0; k<n_atoms; k++ )
+ {
+ if( i == k ) continue;
+
+ rv[0] = pos[0] - atom_pos[k*3 +0];
+ rv[1] = pos[1] - atom_pos[k*3 +1];
+ rv[2] = pos[2] - atom_pos[k*3 +2];
+ r2 = rv[0]*rv[0] + rv[1]*rv[1] + rv[2]*rv[2];
+
+ if( r2 < r*r )
+ {
+ good_point = 0;
+ break;
+ }
+ }
+
+ if( good_point == 1 )
+ {
+ x.push_back( pos[0] );
+ y.push_back( pos[1] );
+ z.push_back( pos[2] );
+ }
+ }
+ }
+
+ // add good vertices to array
+ *n_shell_points = (int)x.size();
+ shell_points = (FLOAT_TYPE *)malloc( sizeof(FLOAT_TYPE)*3*(*n_shell_points) );
+ for( i=0; i<(*n_shell_points); i++ )
+ {
+ shell_points[i*3 +0] = x[i];
+ shell_points[i*3 +1] = y[i];
+ shell_points[i*3 +2] = z[i];
+ }
+
+ // final check; no points should be closer than r from any atom.
+ for( i=0; i<n_atoms; i++ )
+ {
+ //printf( "%f, %f, %f\n", atom_pos[i*3 +0], atom_pos[i*3 +1], atom_pos[i*3 +2] );
+ for( j=0; j<(*n_shell_points); j++ )
+ {
+ //printf( "\t %f, %f, %f\n", shell_points[j*3 +0], shell_points[j*3 +1], shell_points[j*3 +2] );
+ rv[0] = atom_pos[i*3 +0] - shell_points[j*3 +0];
+ rv[1] = atom_pos[i*3 +1] - shell_points[j*3 +1];
+ rv[2] = atom_pos[i*3 +2] - shell_points[j*3 +2];
+ r2 = rv[0]*rv[0] + rv[1]*rv[1] + rv[2]*rv[2];
+
+ if( r2 < (r*r)*0.999 )
+ {
+ printf( "Subdivider::%s(): AAARGH!\n", __func__ );
+ printf( "Point %d is %e from atom %d; minimum is specified as %e\n", j, sqrt(r2), i, r );
+ exit( -1 );
+ }
+ }
+ }
+
+ return shell_points;
+}
+
+// pass 2 random numbers on the range 0 -> 1!
+void UniformSpherePoint( FLOAT_TYPE random_number1, FLOAT_TYPE random_number2, FLOAT_TYPE *vec )
+{
+ /*
+ // original version; this had issues
+
+ double u, theta;
+
+ assert( vec != NULL );
+
+ u = (random_number1-0.5) * 2.0; // u \elem [-1,1]
+ u = sqrt( 1 - u*u );
+ theta = random_number2 * 2.0*M_PI; // theta \elem [0,2pi)
+ vec[0] = u * cos( theta );
+ vec[1] = u * sin( theta );
+ vec[2] = u;
+ */
+
+ // new; tested ok! previous version was wrong.
+ double z, phi, theta, ctheta;
+
+ assert( vec != NULL );
+
+ z = (random_number1 - 0.5) * 2.0;
+ phi = random_number2 * 2.0*M_PI;
+ theta = asin( z );
+ ctheta = cos( theta );
+
+ vec[0] = ctheta * cos(phi);
+ vec[1] = ctheta * sin(phi);
+ vec[2] = z;
+
+}
Property changes on: SwiftApps/CMTS/dimers/src/GeometryUtil.cpp
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/GeometryUtil.h
===================================================================
--- SwiftApps/CMTS/dimers/src/GeometryUtil.h (rev 0)
+++ SwiftApps/CMTS/dimers/src/GeometryUtil.h 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,65 @@
+#ifndef GEOMETRYUTIL_INCLUDED
+
+#include "definitions.h"
+
+#include <assert.h>
+#include <math.h>
+
+/*
+ Where we recieve or pass back data, data type is FLOAT_TYPE. However, function-local variables are double precision.
+*/
+
+#define DOT_PRODUCT( v1, v2 ) ( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] )
+
+FLOAT_TYPE * GenerateSphere( int theta_res, int phi_res, FLOAT_TYPE sphere_radius, int *npoints );
+
+void CartesianToSpherical( FLOAT_TYPE *p1, FLOAT_TYPE *p2, FLOAT_TYPE *r, FLOAT_TYPE *theta, FLOAT_TYPE *phi );
+void SphericalToCartesian( FLOAT_TYPE r, FLOAT_TYPE theta, FLOAT_TYPE phi, FLOAT_TYPE *x, FLOAT_TYPE *y, FLOAT_TYPE *z );
+
+void CentreOfMass( int N, FLOAT_TYPE *coords, FLOAT_TYPE *mass, FLOAT_TYPE *COM );
+
+// rotate on point rotates ... about the specified point! Can pass COM for molecular rotations etc.
+// point = NULL rotates about point (0.0, 0.0, 0.0)
+void RotateAboutPoint( int N, FLOAT_TYPE *xyz, FLOAT_TYPE *point, FLOAT_TYPE *axis, FLOAT_TYPE theta );
+
+// pass 2 random numbers on the range 0 -> 1!
+void UniformSpherePoint( FLOAT_TYPE random_number1, FLOAT_TYPE random_number2, FLOAT_TYPE *vec );
+
+
+FLOAT_TYPE l2_norm_difference( int len, FLOAT_TYPE *v1, FLOAT_TYPE *v2 );
+
+#define INIT_TETRAHEDRON 0
+#define INIT_OCTAHEDRON 1
+#define INIT_ICOSAHEDRON 2
+
+#include <vector>
+
+class Subdivider
+{
+ public:
+
+ int n_vertices, n_faces, n_edges;
+ float *vertices;
+ int *faces;
+
+ int edge_walk;
+ int *start, *end, *midpoint;
+
+ // initial generation stuff
+ void InitTetrahedron();
+ void InitOctahedron();
+ void InitIcosahedron();
+ // refinement of structures
+ int SearchMidpoint( int index_start, int index_end );
+ void Subdivide();
+
+ // enclosing mesh production for a set of atoms, where r is the mesh distance from the atoms
+ FLOAT_TYPE * BuildShell( int n_atoms, FLOAT_TYPE *atom_pos, FLOAT_TYPE r, int *n_shell_points );
+
+ Subdivider( int init_as, int n_subdivisions );
+ ~Subdivider();
+};
+
+
+#define GEOMETRYUTIL_INCLUDED
+#endif
Property changes on: SwiftApps/CMTS/dimers/src/GeometryUtil.h
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/Makefile
===================================================================
--- SwiftApps/CMTS/dimers/src/Makefile (rev 0)
+++ SwiftApps/CMTS/dimers/src/Makefile 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,14 @@
+CC=g++
+
+all: insert_molecules restart2data
+
+insert_molecules:
+ $(CC) -Wall -Wextra -pedantic -O3 insert_molecules.cpp StringUtil.cpp GeometryUtil.cpp Ran1.cpp -o insert_molecules
+ cp insert_molecules ../bin
+
+restart2data:
+ $(CC) -O3 restart2data.cpp -o restart2data
+ cp restart2data ../bin
+
+clean:
+ rm -f restart2data insert_molecules
Added: SwiftApps/CMTS/dimers/src/Ran1.cpp
===================================================================
--- SwiftApps/CMTS/dimers/src/Ran1.cpp (rev 0)
+++ SwiftApps/CMTS/dimers/src/Ran1.cpp 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,58 @@
+#include "Ran1.h"
+
+/*
+ The ran1() function from Numerical Recipes in C: The Art of Scientific Computing ( ISBN: 0-521-43108-5 ).
+ Call with idum as negative number to start, then do not alter idum between calls!
+*/
+#define IA 16807
+#define IM 2147483647
+#define AM (1.0/IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+float ran1( long * idum )
+{
+ int j;
+ long k;
+ static long iy = 0;
+ static long iv[NTAB];
+ float temp;
+
+ if( *idum <= 0 || !iy )
+ {
+ /* initialize */
+ if (-(*idum) < 1) *idum=1; /* prevent idum == 0 */
+ else *idum =-(*idum);
+
+ for( j = NTAB+7; j >= 0; j--)
+ {
+ /* Load the shuffle table ( after 8 warm-ups ) */
+ k = (*idum) / IQ;
+ *idum = IA*(*idum-k*IQ)-IR*k;
+ if( *idum < 0 ) *idum += IM;
+ if( j < NTAB ) iv[j] = *idum;
+ }
+ iy = iv[0];
+ }
+
+ k = (*idum) / IQ; /* Start here when not initializing */
+ *idum=IA*(*idum-k*IQ)-IR*k; /* Compute idum = (IA*idum) % IM without overflows by Schrage's method */
+ if( *idum < 0 ) *idum += IM;
+ j = iy / NDIV; /* Will be in the range 0..NTAB-1. */
+ iy = iv[j]; /* Output previously stored value and refill the shuffle table */
+ iv[j] = *idum;
+ if( (temp = AM*iy) > RNMX ) return RNMX; /* Because users don’t expect end point values. */
+ else return temp;
+}
+#undef IA
+#undef IM
+#undef AM
+#undef IQ
+#undef IR
+#undef NTAB
+#undef NDIV
+#undef EPS
+#undef RNMX
Property changes on: SwiftApps/CMTS/dimers/src/Ran1.cpp
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/Ran1.h
===================================================================
--- SwiftApps/CMTS/dimers/src/Ran1.h (rev 0)
+++ SwiftApps/CMTS/dimers/src/Ran1.h 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,6 @@
+#ifndef RAN1_DEFINED
+
+float ran1( long * idum );
+
+#define RAN1_DEFINED
+#endif
Property changes on: SwiftApps/CMTS/dimers/src/Ran1.h
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/StringUtil.cpp
===================================================================
--- SwiftApps/CMTS/dimers/src/StringUtil.cpp (rev 0)
+++ SwiftApps/CMTS/dimers/src/StringUtil.cpp 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,191 @@
+#include "StringUtil.h"
+
+#define STRINGERROR( msg ) { printf( "ERROR: %s(): %s line %d: \"%s\"\n", __func__, __FILE__, __LINE__, msg ); exit(-1); }
+
+/*
+ Convert null terminated character array into an integer. Checks for bad input string, and produces error where failed.
+*/
+int StringToInt( const char * str, int base, const char * file, int lineno )
+{
+ char *endptr;
+ char message[1024];
+ int returned;
+
+ #ifdef DEBUG
+ if( str == NULL ) STRINGERROR( "NULL string passed as input" );
+ if( file == NULL ) STRINGERROR( "NULL file name passed as input" );
+ #endif
+
+ returned = strtol( str, &endptr, 10 );
+ if( *endptr != '\0' )
+ {
+ sprintf( message, "Unable to convert '%s' (%s, line %d) into an integer of base %d", str, file, lineno, base );
+ STRINGERROR( message );
+ }
+ return returned;
+}
+/*
+ Convert null terminated character array into a double precision floating point. Checks for bad input string, and produces error where failed.
+*/
+double StringToDouble( const char * str, const char * file, int lineno )
+{
+ char *endptr;
+ char message[1024];
+ double returned;
+
+ #ifdef DEBUG
+ if( str == NULL ) STRINGERROR( "NULL string passed" );
+ if( file == NULL ) STRINGERROR( "NULL file name passed" );
+ #endif
+
+ returned = strtod( str, &endptr );
+ if( *endptr != '\0' )
+ {
+ sprintf( message, "Unable to convert '%s' (%s, line %d) into a double", str, file, lineno );
+ STRINGERROR( message );
+ }
+ return returned;
+}
+/*
+ Check if character 'test' is a member of "delimiters"
+*/
+int IsDelim( char test, const char * delimiters, int ndelim )
+{
+ int i;
+
+ #ifdef DEBUG
+ if( delimiters == NULL ) STRINGERROR( "NULL delimiter string passed" );
+ #endif
+
+ for( i=0; i<ndelim; i++ )
+ if( test == delimiters[i] ) return 1;
+ return 0;
+}
+/*
+ Strip whitespace from start and end of string, with user-defined whitespace.
+ Where whitespace pointer is NULL, uses default whitespace of space, tab and newline characters.
+*/
+int StripString( char * str, const char *whitespace )
+{
+ int str_length, whitespace_length, i, j, start, end;
+ const char *default_whitespace = " \t\n";
+ const char *whitespace_ptr;
+
+ #ifdef DEBUG
+ if( str == NULL ) STRINGERROR( "NULL delimiter string passed" );
+ #endif
+
+ if( whitespace != NULL ) whitespace_ptr = whitespace;
+ else whitespace_ptr = default_whitespace;
+ str_length = strlen( str );
+ whitespace_length = strlen( whitespace_ptr );
+
+ /* get index of first non-whitespace character in source string */
+ start = -1;
+ for( i=0; i<str_length && start == -1; i++ )
+ {
+ if( IsDelim( str[i], whitespace_ptr, whitespace_length ) == 0 ) start = i;
+ }
+
+ /* get index of last non-whitespace character in source string */
+ end = -1;
+ for( i=str_length-1; i>=0 && end == -1; i-- )
+ {
+ if( IsDelim( str[i], whitespace_ptr, whitespace_length ) == 0 ) end = i;
+ }
+
+ /* note that we could have an empty string, so check for that. */
+ j = 0;
+ if( start == -1 || end == -1 || end < start )
+ {
+ }
+ else
+ {
+ for( i=start; i<=end; i++ )
+ {
+ str[j] = str[i];
+ j++;
+ }
+ }
+ str[j] = '\0';
+
+ return strlen( str );
+}
+/*
+ Chops up "str", and stores the substrings in the array "pointers". If too many substrings are present to store in
+ the "pointers" array, this routine will just do as many as possible.
+*/
+int TokenizeString( char * str, const char * delimiters, char ** pointers, int maxptrs )
+{
+ int ntoks, len, ndelim, i, j;
+
+ #ifdef DEBUG
+ if( str == NULL ) STRINGERROR( "NULL string passed" );
+ if( delimiters == NULL ) STRINGERROR( "NULL delimiters string passed" );
+ if( pointers == NULL ) STRINGERROR( "NULL pointers array passed" );
+ if( maxptrs < 1 ) STRINGERROR( "maxptrs < 1" );
+ #endif
+
+ len = strlen( str );
+ ndelim = strlen( delimiters );
+ ntoks = 0;
+
+ j = 0; /* index into current token. */
+ for( i=0; i<len; i++ )
+ {
+ if( IsDelim( str[i], delimiters, ndelim ) == 1 )
+ {
+ /* If we've been generating a valid token and hit whitespace, increase ntokens and reset token character index */
+ if( j > 0 ) ntoks++;
+ j = 0;
+ if( ntoks >= maxptrs ) break;
+ }
+ else
+ {
+ pointers[ntoks][j] = str[i];
+ pointers[ntoks][j+1] = '\0';
+ j++;
+ }
+ }
+ /* We may have a trailing token at the end of the string */
+ if( j > 0 ) ntoks++;
+ return ntoks;
+}
+
+int TokenizeString( const char * source, const char * delimiters, std::vector< std::string > &results )
+{
+ int ntoks;
+ char buffer[1024], tokens[100][1024], *token_ptrs[100];
+ std::string str;
+
+ for( int i=0; i<100; i++ ) token_ptrs[i] = tokens[i];
+
+ strcpy( buffer, source );
+
+ results.clear();
+
+ ntoks = TokenizeString( buffer, delimiters, token_ptrs, 100 );
+ for( int i=0; i<ntoks; i++ )
+ {
+ str = tokens[i];
+ results.push_back( str );
+ }
+
+ return ntoks;
+}
+int TokenizeLine( FILE * f, const char *delimiters, char *line_buffer, std::vector<std::string> &results )
+{
+ #ifdef DEBUG
+ if( f == NULL ) STRINGERROR( "NULL FILE* passed" );
+ if( delimiters == NULL ) STRINGERROR( "NULL delimiters string passed" );
+ if( line_buffer == NULL ) STRINGERROR( "NULL line_buffer pointer passed" );
+ #endif
+
+ line_buffer[0] = '\0';
+ if( fgets( line_buffer, 1024, f) == NULL ) return -1;
+ return TokenizeString( line_buffer, delimiters, results );
+}
+
+
+#undef STRINGERROR
+
Property changes on: SwiftApps/CMTS/dimers/src/StringUtil.cpp
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/StringUtil.h
===================================================================
--- SwiftApps/CMTS/dimers/src/StringUtil.h (rev 0)
+++ SwiftApps/CMTS/dimers/src/StringUtil.h 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,19 @@
+#ifndef STRING_UTIL_DEFINED
+
+#include <string>
+#include <vector>
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+int StringToInt( const char * str, int base, const char * file, int lineno );
+double StringToDouble( const char * str, const char * file, int lineno );
+int IsDelim( char test, const char * delimiters, int ndelim );
+int StripString( char * str, const char *whitespace );
+int TokenizeString( char * str, const char * delimiters, char ** pointers, int maxptrs );
+int TokenizeString( const char * source, const char * delimiters, std::vector< std::string > &results );
+int TokenizeLine( FILE * f, const char *delimiters, char *line_buffer, std::vector<std::string> &results );
+
+#define STRING_UTIL_DEFINED
+#endif
Property changes on: SwiftApps/CMTS/dimers/src/StringUtil.h
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/definitions.h
===================================================================
--- SwiftApps/CMTS/dimers/src/definitions.h (rev 0)
+++ SwiftApps/CMTS/dimers/src/definitions.h 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1 @@
+#define FLOAT_TYPE float
Property changes on: SwiftApps/CMTS/dimers/src/definitions.h
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/insert_molecules
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/src/insert_molecules
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/src/insert_molecules.cpp
===================================================================
--- SwiftApps/CMTS/dimers/src/insert_molecules.cpp (rev 0)
+++ SwiftApps/CMTS/dimers/src/insert_molecules.cpp 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,700 @@
+#include <stdio.h>
+#include <math.h>
+
+#include <map>
+#include <vector>
+#include <string>
+
+#include "StringUtil.h"
+#include "GeometryUtil.h"
+#include "Ran1.h"
+
+/*
+ This inserts one or more copies of a lammps config file into another lammps config file.
+ Bonding info and molecular assignments etc are preserved, and incremented as appropriate.
+
+ The program assumes you're inserting into a sphere, and will insert in the whole sphere
+ volume or onto the surface as specified. Also specified is a cutoff distance to any existing
+ site in the simulation to avoid steric clashes.
+
+ g++ -Wall -Wextra -pedantic -O3 config_insert.cpp StringUtil.cpp GeometryUtil.cpp Ran1.cpp -o config_insert
+*/
+
+struct atom_data
+{
+ int serial, mol, type;
+ double x, y, z;
+ double vx, vy, vz;
+};
+struct bond_data
+{
+ int serial, type, i, j;
+};
+struct improper_data
+{
+ int serial, type, i, j, k, l;
+};
+struct sim_data
+{
+ double minx,maxx;
+ double miny,maxy;
+ double minz,maxz;
+
+ int n_atoms, n_bonds, n_impropers;
+ int n_atom_types, n_bond_types, n_improper_types;
+
+ std::map<int, double> masses;
+ std::map<int, struct atom_data> atoms;
+ std::map<int, struct bond_data> bonds;
+ std::map<int, struct improper_data> impropers;
+};
+void add_atom( struct sim_data &sd, int type, int mol, double x, double y, double z, double vx, double vy, double vz )
+{
+ struct atom_data ad;
+
+ if( sd.masses.find( type ) == sd.masses.end() )
+ {
+ printf( "%s(): unable to find atom type %d in defined masses.\n", __func__, type );
+ return;
+ }
+
+ sd.n_atoms++;
+
+ ad.serial = sd.n_atoms;
+ ad.mol = mol;
+ ad.type = type;
+ ad.x = x;
+ ad.y = y;
+ ad.z = z;
+ ad.vx = vx;
+ ad.vy = vy;
+ ad.vz = vz;
+
+ sd.atoms[ ad.serial ] = ad;
+}
+void add_bond( struct sim_data &sd, int type, int i, int j )
+{
+ struct bond_data bd;
+
+ if( sd.atoms.find( i ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom index %d in defined atoms.\n", __func__, i );
+ return;
+ }
+ if( sd.atoms.find( j ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom index %d in defined atoms.\n", __func__, j );
+ return;
+ }
+
+ sd.n_bonds++;
+
+ bd.serial = sd.n_bonds;
+ bd.type = type;
+ bd.i = i;
+ bd.j = j;
+
+ sd.bonds[ bd.serial ] = bd;
+}
+void add_improper( struct sim_data &sd, int type, int i, int j, int k, int l )
+{
+ struct improper_data id;
+
+ sd.n_impropers++;
+
+ id.serial = sd.n_impropers;
+ id.type = type;
+ id.i = i;
+ id.j = j;
+ id.k = k;
+ id.l = l;
+
+ sd.impropers[ id.serial ] = id;
+}
+int get_largest_molecule_id( struct sim_data &sd )
+{
+ std::map<int, struct atom_data>::iterator it;
+ int largest;
+
+ largest = 0;
+
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ if( it->second.mol > largest ) largest = it->second.mol;
+ }
+ return largest;
+}
+
+void load_data(
+ FILE * f,
+ struct sim_data &sd
+ )
+{
+ char line_buffer[1024];
+ const char *delimiters = " \t\n";
+ int line_no;
+
+ std::string state;
+ std::vector<std::string> tokens;
+ std::map<std::string,int> states;
+
+ std::map<int, struct atom_data>::iterator a_it;
+ std::map<int, struct bond_data>::iterator b_it;
+ std::map<int, struct improper_data>::iterator i_it;
+ std::map<int,int> int_map;
+
+ states["Masses"] = 1;
+ states["Atoms"] = 1;
+ states["Velocities"] = 1;
+ states["Bonds"] = 1;
+ states["Impropers"] = 1;
+ // following are in case the coeffs are included in the file, so we don't get confused, eg "Bond Coeffs"
+ states["Bond"] = 1;
+ states["Improper"] = 1;
+
+ state = "";
+
+ sd.minx = sd.maxx = 0.0;
+ sd.miny = sd.maxy = 0.0;
+ sd.minz = sd.maxz = 0.0;
+ sd.n_atoms = sd.n_atom_types = 0;
+ sd.n_bonds = sd.n_bond_types = 0;
+ sd.n_impropers = sd.n_improper_types = 0;
+ sd.masses.clear();
+ sd.atoms.clear();
+ sd.bonds.clear();
+ sd.impropers.clear();
+
+ TokenizeLine( f, delimiters, line_buffer, tokens ); // skip leading line.
+
+ while( TokenizeLine( f, delimiters, line_buffer, tokens ) != -1 )
+ {
+ line_no++;
+
+ if( tokens.size() < 1 ) continue;
+
+ if( states.find( tokens[0] ) != states.end() )
+ {
+ //printf( "State switch: '%s' => '%s'\n", state.c_str(), tokens[0].c_str() );
+ state = tokens[0];
+ continue;
+ }
+
+ if( state == "" )
+ {
+ // load cell etc, put into metadata
+ if( tokens.size() > 3 )
+ {
+ if( tokens[2] == "xlo" && tokens[3] == "xhi" )
+ {
+ sd.minx = atof( tokens[0].c_str() );
+ sd.maxx = atof( tokens[1].c_str() );
+ }
+ if( tokens[2] == "ylo" && tokens[3] == "yhi" )
+ {
+ sd.miny = atof( tokens[0].c_str() );
+ sd.maxy = atof( tokens[1].c_str() );
+ }
+ if( tokens[2] == "zlo" && tokens[3] == "zhi" )
+ {
+ sd.minz = atof( tokens[0].c_str() );
+ sd.maxz = atof( tokens[1].c_str() );
+ }
+ }
+ else if( tokens.size() == 2 && tokens[1] == "atoms" ) { sd.n_atoms = atoi( tokens[0].c_str() ); }
+ else if( tokens.size() == 2 && tokens[1] == "bonds" ) { sd.n_bonds = atoi( tokens[0].c_str() ); }
+ else if( tokens.size() == 2 && tokens[1] == "impropers" ) { sd.n_impropers = atoi( tokens[0].c_str() ); }
+ else if( tokens.size() == 3 && tokens[1] == "atom" ) { sd.n_atom_types = atoi( tokens[0].c_str() ); }
+ else if( tokens.size() == 3 && tokens[1] == "bond" ) { sd.n_bond_types = atoi( tokens[0].c_str() ); }
+ else if( tokens.size() == 3 && tokens[1] == "improper" ) { sd.n_improper_types = atoi( tokens[0].c_str() ); }
+ }
+ else if( state == "Masses" )
+ {
+ sd.masses[ atoi( tokens[0].c_str() ) ] = atof( tokens[1].c_str() );
+ }
+ else if( state == "Atoms" )
+ {
+ struct atom_data ad;
+
+ ad.serial = atoi( tokens[0].c_str() );
+ ad.mol = atoi( tokens[1].c_str() );
+ ad.type = atoi( tokens[2].c_str() );
+
+ ad.x = atof( tokens[3].c_str() );
+ ad.y = atof( tokens[4].c_str() );
+ ad.z = atof( tokens[5].c_str() );
+
+ ad.vx = 0.0;
+ ad.vy = 0.0;
+ ad.vz = 0.0;
+
+ // assume zero velocoties by default ; will be overwritten in Velocities section, where exists
+
+ sd.atoms[ ad.serial ] = ad;
+ }
+ else if( state == "Velocities" )
+ {
+ int serial;
+
+ serial = atoi( tokens[0].c_str() );
+ if( sd.atoms.find( serial ) == sd.atoms.end() )
+ {
+ printf( "%s(): '%s' : unable to find atom with serial %d\n", __func__, state.c_str(), serial );
+ exit( -1 );
+ }
+
+ sd.atoms[serial].vx = atof( tokens[1].c_str() );
+ sd.atoms[serial].vy = atof( tokens[2].c_str() );
+ sd.atoms[serial].vz = atof( tokens[3].c_str() );
+ }
+ else if( state == "Bonds" )
+ {
+ struct bond_data bd;
+
+ bd.serial = atoi( tokens[0].c_str() );
+ bd.type = atoi( tokens[1].c_str() );
+ bd.i = atoi( tokens[2].c_str() );
+ bd.j = atoi( tokens[3].c_str() );
+
+ sd.bonds[ bd.serial ] = bd;
+
+ //printf( "[%d] Bond: %d, %d, %d, %d\n", (int)sd.bonds.size(), bd.serial, bd.type, bd.i, bd.j );
+ }
+ else if( state == "Impropers" )
+ {
+ struct improper_data id;
+
+ id.serial = atoi( tokens[0].c_str() );
+ id.type = atoi( tokens[1].c_str() );
+ id.i = atoi( tokens[2].c_str() );
+ id.j = atoi( tokens[3].c_str() );
+ id.k = atoi( tokens[4].c_str() );
+ id.l = atoi( tokens[5].c_str() );
+
+ sd.impropers[ id.serial ] = id;
+ //printf( "[%d] Improper: %d, %d, %d, %d, %d, %d\n", (int)sd.impropers.size(), id.serial, id.type, id.i, id.j, id.k, id.l );
+ }
+ }
+
+ // error tests, and build system info.
+
+ for( a_it = sd.atoms.begin(); a_it != sd.atoms.end(); a_it++ )
+ {
+ if( sd.masses.find( a_it->second.type ) == sd.masses.end() )
+ {
+ printf( "%s(): unable to find atom mass for type %d\n", __func__, a_it->second.type );
+ exit( -1 );
+ }
+ }
+ if( sd.n_atoms != (int)sd.atoms.size() )
+ {
+ printf( "%s(): system declares %d atoms, but only have data for %d\n", __func__, sd.n_atoms, (int)sd.atoms.size() );
+ exit( -1 );
+ }
+
+ for( b_it = sd.bonds.begin(); b_it != sd.bonds.end(); b_it++ )
+ {
+ if( sd.atoms.find( b_it->second.i ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in bond definition\n", __func__, b_it->second.i );
+ exit( -1 );
+ }
+ if( sd.atoms.find( b_it->second.j ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in bond definition\n", __func__, b_it->second.j );
+ exit( -1 );
+ }
+ if( b_it->second.type > sd.n_bond_types )
+ {
+ printf( "%s(): bond type %d is larger than declared number of types %d\n", __func__, b_it->second.type, sd.n_bond_types );
+ exit( -1 );
+ }
+ //printf( "Bond: %d %d %d %d\n", b_it->second.serial, b_it->second.type, b_it->second.i, b_it->second.j );
+ }
+ if( sd.n_bonds != (int)sd.bonds.size() )
+ {
+ printf( "%s(): system declares %d bonds, but only have data for %d\n", __func__, sd.n_bonds, (int)sd.bonds.size() );
+ exit( -1 );
+ }
+
+ for( i_it = sd.impropers.begin(); i_it != sd.impropers.end(); i_it++ )
+ {
+ if( sd.atoms.find( i_it->second.i ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in improper definition\n", __func__, i_it->second.i );
+ exit( -1 );
+ }
+ if( sd.atoms.find( i_it->second.j ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in improper definition\n", __func__, i_it->second.j );
+ exit( -1 );
+ }
+ if( sd.atoms.find( i_it->second.k ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in improper definition\n", __func__, i_it->second.k );
+ exit( -1 );
+ }
+ if( sd.atoms.find( i_it->second.l ) == sd.atoms.end() )
+ {
+ printf( "%s(): unable to find atom serial %d in improper definition\n", __func__, i_it->second.l );
+ exit( -1 );
+ }
+ if( i_it->second.type > sd.n_improper_types )
+ {
+ printf( "%s(): improper type %d is larger than declared number of types %d\n", __func__, i_it->second.type, sd.n_improper_types );
+ exit( -1 );
+ }
+ }
+ if( sd.n_impropers != (int)sd.impropers.size() )
+ {
+ printf( "%s(): system declares %d impropers, but only have data for %d\n", __func__, sd.n_impropers, (int)sd.impropers.size() );
+ exit( -1 );
+ }
+
+}
+
+void save_data(
+ FILE * f,
+ struct sim_data &sd
+ )
+{
+ std::map<int, double>::iterator m_it;
+ std::map<int, struct atom_data>::iterator a_it;
+ std::map<int, struct bond_data>::iterator b_it;
+ std::map<int, struct improper_data>::iterator i_it;
+
+ fprintf( f, "Autogenerated!\n" );
+ fprintf( f, "\n" );
+ fprintf( f, "%d atoms\n", sd.n_atoms );
+ fprintf( f, "%d bonds\n", sd.n_bonds );
+ fprintf( f, "%d impropers\n", sd.n_impropers );
+ fprintf( f, "\n" );
+ fprintf( f, "%d atom types\n", sd.n_atom_types );
+ fprintf( f, "%d bond types\n", sd.n_bond_types );
+ fprintf( f, "%d improper types\n", sd.n_improper_types );
+ fprintf( f, "\n" );
+ fprintf( f, "%+g %+g xlo xhi\n", sd.minx, sd.maxx );
+ fprintf( f, "%+g %+g ylo yhi\n", sd.miny, sd.maxy );
+ fprintf( f, "%+g %+g zlo zhi\n", sd.minz, sd.maxz );
+ fprintf( f, "\n" );
+ fprintf( f, "Masses\n" );
+ fprintf( f, "\n" );
+ for( m_it = sd.masses.begin(); m_it != sd.masses.end(); m_it++ ) fprintf( f, "%d %g\n", m_it->first, m_it->second );
+ fprintf( f, "\n" );
+ fprintf( f, "Atoms\n" );
+ fprintf( f, "\n" );
+ for( a_it = sd.atoms.begin(); a_it != sd.atoms.end(); a_it++ ) fprintf( f, "%d %d %d %g %g %g\n", a_it->first, a_it->second.mol, a_it->second.type, a_it->second.x, a_it->second.y, a_it->second.z );
+ fprintf( f, "\n" );
+ fprintf( f, "Velocities\n" );
+ fprintf( f, "\n" );
+ for( a_it = sd.atoms.begin(); a_it != sd.atoms.end(); a_it++ ) fprintf( f, "%d %g %g %g\n", a_it->first, a_it->second.vx, a_it->second.vy, a_it->second.vz );
+ fprintf( f, "\n" );
+ fprintf( f, "Bonds\n" );
+ fprintf( f, "\n" );
+ for( b_it = sd.bonds.begin(); b_it != sd.bonds.end(); b_it++ ) fprintf( f, "%d %d %d %d\n", b_it->first, b_it->second.type, b_it->second.i, b_it->second.j );
+ fprintf( f, "\n" );
+ fprintf( f, "Impropers\n" );
+ fprintf( f, "\n" );
+ for( i_it = sd.impropers.begin(); i_it != sd.impropers.end(); i_it++ ) fprintf( f, "%d %d %d %d %d %d\n", i_it->first, i_it->second.type, i_it->second.i, i_it->second.j, i_it->second.k, i_it->second.l );
+ fprintf( f, "\n" );
+}
+
+double reset_cog( struct sim_data &sd )
+{
+ std::map<int, struct atom_data>::iterator it;
+ double max_extent_from_cog_sq;
+ double COG[3], COG2[3];
+
+ COG[0] = COG[1] = COG[2] = 0.0;
+ COG2[0] = COG2[1] = COG2[2] = 0.0;
+
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ COG[0] += it->second.x;
+ COG[1] += it->second.y;
+ COG[2] += it->second.z;
+ }
+ COG[0] /= sd.n_atoms;
+ COG[1] /= sd.n_atoms;
+ COG[2] /= sd.n_atoms;
+ printf( "COM was: %g, %g, %g\n", COG[0], COG[1], COG[2] );
+ max_extent_from_cog_sq = -1.0;
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ double dr[3], extent_sq;
+
+ it->second.x -= COG[0];
+ it->second.y -= COG[1];
+ it->second.z -= COG[2];
+
+ COG2[0] += it->second.x;
+ COG2[1] += it->second.y;
+ COG2[2] += it->second.z;
+
+ // COG now at origin, or should be, so just use the coords.
+ dr[0] = it->second.x;
+ dr[1] = it->second.y;
+ dr[2] = it->second.z;
+ extent_sq = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
+ if( extent_sq > max_extent_from_cog_sq )
+ {
+ max_extent_from_cog_sq = extent_sq;
+ }
+ }
+ COG2[0] /= sd.n_atoms;
+ COG2[1] /= sd.n_atoms;
+ COG2[2] /= sd.n_atoms;
+ printf( "COM now: %g, %g, %g\n", COG2[0], COG2[1], COG2[2] );
+
+ return sqrt( max_extent_from_cog_sq );
+}
+
+void add_fragment( struct sim_data &sd, struct sim_data &frag )
+{
+ int start_mol_id, start_atom_serial;
+
+ std::map<int, struct atom_data>::iterator a_it;
+ std::map<int, struct bond_data>::iterator b_it;
+ std::map<int, struct improper_data>::iterator i_it;
+
+ start_mol_id = get_largest_molecule_id( sd );
+ start_atom_serial = sd.n_atoms;
+
+ for( a_it = frag.atoms.begin(); a_it != frag.atoms.end(); a_it++ )
+ {
+ add_atom( sd,
+ a_it->second.type,
+ a_it->second.mol+start_mol_id,
+ a_it->second.x,
+ a_it->second.y,
+ a_it->second.z,
+ a_it->second.vx,
+ a_it->second.vy,
+ a_it->second.vz );
+ }
+ for( b_it = frag.bonds.begin(); b_it != frag.bonds.end(); b_it++ )
+ {
+ add_bond( sd, b_it->second.type, b_it->second.i+start_atom_serial, b_it->second.j+start_atom_serial );
+ }
+ for( i_it = frag.impropers.begin(); i_it != frag.impropers.end(); i_it++ )
+ {
+ add_improper( sd, i_it->second.type, i_it->second.i+start_atom_serial, i_it->second.j+start_atom_serial, i_it->second.k+start_atom_serial, i_it->second.l+start_atom_serial );
+ }
+}
+
+int check_collisions( struct sim_data &sd, struct sim_data &frag, float min_separation )
+{
+ std::map<int, struct atom_data>::iterator i_it;
+ std::map<int, struct atom_data>::iterator j_it;
+
+ for( i_it = frag.atoms.begin(); i_it != frag.atoms.end(); i_it++ )
+ {
+ for( j_it = sd.atoms.begin(); j_it != sd.atoms.end(); j_it++ )
+ {
+ double dr[3], r_sq;
+
+ dr[0] = i_it->second.x - j_it->second.x;
+ dr[1] = i_it->second.y - j_it->second.y;
+ dr[2] = i_it->second.z - j_it->second.z;
+ r_sq = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
+
+ if( r_sq < min_separation*min_separation ) return 1;
+ }
+ }
+
+ return 0;
+}
+
+void sim_coords_to_xyz( struct sim_data &sd, FLOAT_TYPE *xyz )
+{
+ std::map<int, struct atom_data>::iterator it;
+ int i;
+
+ i = 0;
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ xyz[i*3 +0] = it->second.x;
+ xyz[i*3 +1] = it->second.y;
+ xyz[i*3 +2] = it->second.z;
+ i++;
+ }
+}
+void xyz_to_sim_coords( FLOAT_TYPE *xyz, struct sim_data &sd )
+{
+ std::map<int, struct atom_data>::iterator it;
+ int i;
+
+ i = 0;
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ it->second.x = xyz[i*3+0];
+ it->second.y = xyz[i*3+1];
+ it->second.z = xyz[i*3+2];
+ i++;
+ }
+}
+void save_xyz( const char *fpath, int N, FLOAT_TYPE *xyz )
+{
+ FILE * f;
+
+ f = fopen( fpath, "w" );
+ fprintf( f, "%d\n", N );
+ fprintf( f, "Comment line\n" );
+ for( int i=0; i<N; i++ )
+ {
+ fprintf( f, "P %g %g %g\n", xyz[i*3+0], xyz[i*3+1], xyz[i*3+2] );
+ }
+ fclose( f );
+}
+void save_sim_as_xyz( const char *fpath, struct sim_data &sd )
+{
+ std::map<int, struct atom_data>::iterator it;
+ FILE * f;
+
+ f = fopen( fpath, "w" );
+ fprintf( f, "%d\n", (int)sd.atoms.size() );
+ fprintf( f, "Comment line\n" );
+
+ for( it = sd.atoms.begin(); it != sd.atoms.end(); it++ )
+ {
+ fprintf( f, "P %g %g %g\n", it->second.x, it->second.y, it->second.z );
+ }
+
+ fclose( f );
+}
+
+int main( int argc, char **argv )
+{
+ struct sim_data sim, fragment, frag;
+ int ncopies, insert_on_surface;
+ double radius, min_sep, max_extent_from_cog;
+ long ran1_seed;
+ FLOAT_TYPE *original_xyz, *xyz;
+
+ FILE *f;
+ const char *sim_fpath, *fragment_fpath, *output_fpath;
+
+ if( argc < 8 )
+ {
+ printf( "Usage: %s <sim data> <fragment data> <ncopies of fragment> <sphere radius> <insert_on_surface> <min_separation> <output file>\n", argv[0] );
+ exit( -1 );
+ }
+
+ sim_fpath = argv[1];
+ fragment_fpath = argv[2];
+ ncopies = atoi( argv[3] );
+ radius = atof( argv[4] );
+ insert_on_surface = atoi( argv[5] );
+ min_sep = atof( argv[6] );
+ output_fpath = argv[7];
+
+ // load target sim
+ if( (f = fopen( sim_fpath, "r" )) == NULL )
+ {
+ printf( "Unable to open file '%s'\n", sim_fpath );
+ exit( EXIT_FAILURE );
+ }
+ load_data( f, sim );
+ fclose( f );
+
+ // load fragment to add to target sim
+ if( (f = fopen( fragment_fpath, "r" )) == NULL )
+ {
+ printf( "Unable to open file '%s'\n", fragment_fpath );
+ exit( -1 );
+ }
+ load_data( f, fragment );
+ fclose( f );
+
+ max_extent_from_cog = reset_cog( fragment );
+
+ printf( "Target sim has %d atoms (%d molecules)\n", sim.n_atoms, get_largest_molecule_id( sim ) );
+ printf( "Fragment sim has %d atoms (%d molecules), max extent from COG is %g\n", fragment.n_atoms, get_largest_molecule_id( fragment ), max_extent_from_cog );
+ printf( "%d copies of fragment data to insert.\n", ncopies );
+ printf( "sphere radius is %g.\n", radius );
+
+ original_xyz = new FLOAT_TYPE[ fragment.atoms.size()*3 ];
+ xyz = new FLOAT_TYPE[ fragment.atoms.size()*3 ];
+
+ sim_coords_to_xyz( fragment, original_xyz );
+ frag = fragment;
+
+ for( int i=0; i<ncopies; i++ )
+ {
+ int collided;
+
+ if( ncopies >= 10 && i % (ncopies/10) == 0 ) { printf( "Adding %d...\n", i+1 ); save_sim_as_xyz( "temp.xyz", sim ); }
+
+ collided = 1;
+ for( int j=0; j<1e6; j++ )
+ {
+ FLOAT_TYPE origin[3], COG[3], axis[3], theta;
+ FLOAT_TYPE r1, r2, distance_from_origin;
+
+ memcpy( xyz, original_xyz, sizeof(FLOAT_TYPE)*frag.atoms.size()*3 );
+
+ // generate random axis of rotation
+ r1 = ran1( &ran1_seed );
+ r2 = ran1( &ran1_seed );
+ UniformSpherePoint( r1, r2, axis );
+
+ // generate random rotation about this axis
+ theta = (ran1( &ran1_seed )-0.5) * 2.0*M_PI;
+
+ // generate random COG offset onto surface of sphere
+ r1 = ran1( &ran1_seed );
+ r2 = ran1( &ran1_seed );
+ UniformSpherePoint( r1, r2, COG );
+ if( insert_on_surface == 1 )
+ {
+ distance_from_origin = radius - max_extent_from_cog;
+ }
+ else
+ {
+ distance_from_origin = ran1( &ran1_seed ) * (radius - max_extent_from_cog);
+ }
+ COG[0] *= distance_from_origin;
+ COG[1] *= distance_from_origin;
+ COG[2] *= distance_from_origin;
+
+ origin[0] = origin[1] = origin[2] = 0.0;
+
+ RotateAboutPoint( (int)frag.atoms.size(), xyz, origin, axis, theta );
+ for( size_t i=0; i<frag.atoms.size(); i++ )
+ {
+ xyz[i*3+0] += COG[0];
+ xyz[i*3+1] += COG[1];
+ xyz[i*3+2] += COG[2];
+ }
+ xyz_to_sim_coords( xyz, frag );
+
+ // test for collisions
+ if( check_collisions( sim, frag, min_sep ) == 0 )
+ {
+ collided = 0;
+ break;
+ }
+ }
+ if( collided == 1 )
+ {
+ printf( "Unable to insert fragment.\n" );
+ exit( -1 );
+ }
+
+ add_fragment( sim, frag );
+ }
+
+ delete[] original_xyz;
+ delete[] xyz;
+
+ // save target sim
+ if( (f = fopen( output_fpath, "w" )) == NULL )
+ {
+ printf( "Unable to open file '%s'\n", output_fpath );
+ exit( EXIT_FAILURE );
+ }
+ save_data( f, sim );
+ fclose( f );
+
+ exit( EXIT_SUCCESS );
+}
Property changes on: SwiftApps/CMTS/dimers/src/insert_molecules.cpp
___________________________________________________________________
Added: svn:executable
+ *
Added: SwiftApps/CMTS/dimers/src/restart2data
===================================================================
(Binary files differ)
Property changes on: SwiftApps/CMTS/dimers/src/restart2data
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: SwiftApps/CMTS/dimers/src/restart2data.cpp
===================================================================
--- SwiftApps/CMTS/dimers/src/restart2data.cpp (rev 0)
+++ SwiftApps/CMTS/dimers/src/restart2data.cpp 2012-08-08 15:12:55 UTC (rev 5893)
@@ -0,0 +1,3684 @@
+/* -----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ www.cs.sandia.gov/~sjplimp/lammps.html
+ Steve Plimpton, sjplimp at sandia.gov, Sandia National Laboratories
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------ */
+
+// Convert a LAMMPS binary restart file into an ASCII text data file
+//
+// Syntax: restart2data restart-file data-file (input-file)
+// restart-file and data-file are mandatory
+// input-file is optional
+// if specified it will contain LAMMPS input script commands
+// for mass and force field info
+// only a few force field styles support this option
+//
+// this serial code must be compiled on a platform that can read the binary
+// restart file since binary formats are not compatible across all platforms
+// restart-file can have a '%' character to indicate a multiproc restart
+// file as written by LAMMPS
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+#define MAX_GROUP 32
+#define PI (4.0*atan(1.0))
+
+// these should match settings in src/lmptype.h
+
+#include "stdint.h"
+#define __STDC_FORMAT_MACROS
+#include "inttypes.h"
+
+typedef int tagint;
+typedef int64_t bigint;
+#define BIGINT_FORMAT "%" PRId64
+
+// same as write_restart.cpp
+
+enum{VERSION,SMALLINT,TAGINT,BIGINT,
+ UNITS,NTIMESTEP,DIMENSION,NPROCS,PROCGRID_0,PROCGRID_1,PROCGRID_2,
+ NEWTON_PAIR,NEWTON_BOND,XPERIODIC,YPERIODIC,ZPERIODIC,
+ BOUNDARY_00,BOUNDARY_01,BOUNDARY_10,BOUNDARY_11,BOUNDARY_20,BOUNDARY_21,
+ ATOM_STYLE,NATOMS,NTYPES,
+ NBONDS,NBONDTYPES,BOND_PER_ATOM,
+ NANGLES,NANGLETYPES,ANGLE_PER_ATOM,
+ NDIHEDRALS,NDIHEDRALTYPES,DIHEDRAL_PER_ATOM,
+ NIMPROPERS,NIMPROPERTYPES,IMPROPER_PER_ATOM,
+ BOXLO_0,BOXHI_0,BOXLO_1,BOXHI_1,BOXLO_2,BOXHI_2,
+ SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3,
+ SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3,
+ XY,XZ,YZ};
+enum{MASS};
+enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
+
+static const char * const cg_type_list[] =
+ {"none", "lj9_6", "lj12_4", "lj12_6"};
+
+// ---------------------------------------------------------------------
+// Data class to hold problem
+// ---------------------------------------------------------------------
+
+class Data {
+ public:
+
+ // global settings
+
+ char *version;
+ int size_smallint,size_tagint,size_bigint;
+ bigint ntimestep;
+ int nprocs;
+ char *unit_style;
+ int dimension;
+ int px,py,pz;
+ int newton_pair,newton_bond;
+ int xperiodic,yperiodic,zperiodic;
+ int boundary[3][2];
+
+ char *atom_style;
+ int style_angle,style_atomic,style_bond,style_charge,style_dipole;
+ int style_ellipsoid,style_full;
+ int style_hybrid,style_molecular,style_peri,style_sphere;
+
+ bigint natoms;
+ bigint nellipsoids;
+ bigint nbonds,nangles,ndihedrals,nimpropers;
+ int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
+ int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
+ int triclinic;
+
+ double xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz;
+ double special_lj[4],special_coul[4];
+
+ double cut_lj_global,cut_coul_global,kappa;
+ int offset_flag,mix_flag;
+
+ // force fields
+
+ char *pair_style,*bond_style,*angle_style,*dihedral_style,*improper_style;
+
+ double *pair_born_A,*pair_born_rho,*pair_born_sigma;
+ double *pair_born_C,*pair_born_D;
+ double *pair_buck_A,*pair_buck_rho,*pair_buck_C;
+ double *pair_colloid_A12,*pair_colloid_sigma;
+ double *pair_colloid_d1,*pair_colloid_d2;
+ double *pair_dipole_epsilon,*pair_dipole_sigma;
+ double *pair_dpd_a0,*pair_dpd_gamma;
+ double *pair_charmm_epsilon,*pair_charmm_sigma;
+ double *pair_charmm_eps14,*pair_charmm_sigma14;
+ double *pair_class2_epsilon,*pair_class2_sigma;
+ double *pair_gb_epsilon,*pair_gb_sigma;
+ double *pair_gb_epsa,*pair_gb_epsb,*pair_gb_epsc;
+ double *pair_lj_epsilon,*pair_lj_sigma;
+ double **pair_cg_epsilon,**pair_cg_sigma;
+ int **pair_cg_cmm_type, **pair_setflag;
+ double **pair_cut_coul, **pair_cut_lj;
+ double *pair_ljexpand_epsilon,*pair_ljexpand_sigma,*pair_ljexpand_shift;
+ double *pair_ljgromacs_epsilon,*pair_ljgromacs_sigma;
+ double *pair_ljsmooth_epsilon,*pair_ljsmooth_sigma;
+ double *pair_morse_d0,*pair_morse_alpha,*pair_morse_r0;
+ double *pair_soft_A;
+ double *pair_yukawa_A;
+
+ double *bond_class2_r0,*bond_class2_k2,*bond_class2_k3,*bond_class2_k4;
+ double *bond_fene_k,*bond_fene_r0,*bond_fene_epsilon,*bond_fene_sigma;
+ double *bond_feneexpand_k,*bond_feneexpand_r0;
+ double *bond_feneexpand_epsilon,*bond_feneexpand_sigma;
+ double *bond_feneexpand_shift;
+ double *bond_harmonic_k,*bond_harmonic_r0;
+ double *bond_harmonicshift_umin,*bond_harmonicshift_r0,
+ *bond_harmonicshift_rc;
+ double *bond_harmonicshiftcut_umin,*bond_harmonicshiftcut_r0,
+ *bond_harmonicshiftcut_rc;
+ double *bond_morse_d0,*bond_morse_alpha,*bond_morse_r0;
+ double *bond_nonlinear_epsilon,*bond_nonlinear_r0,*bond_nonlinear_lamda;
+ double *bond_quartic_k,*bond_quartic_b1,*bond_quartic_b2;
+ double *bond_quartic_rc,*bond_quartic_u0;
+
+ double *angle_charmm_k,*angle_charmm_theta0;
+ double *angle_charmm_k_ub,*angle_charmm_r_ub;
+ double *angle_class2_theta0;
+ double *angle_class2_k2,*angle_class2_k3,*angle_class2_k4;
+ double *angle_class2_bb_k,*angle_class2_bb_r1,*angle_class2_bb_r2;
+ double *angle_class2_ba_k1,*angle_class2_ba_k2;
+ double *angle_class2_ba_r1,*angle_class2_ba_r2;
+ double *angle_cosine_k;
+ double *angle_cosine_squared_k,*angle_cosine_squared_theta0;
+ double *angle_harmonic_k,*angle_harmonic_theta0;
+ double *angle_cg_cmm_epsilon,*angle_cg_cmm_sigma;
+ double *angle_cosineshift_umin,*angle_cosineshift_sint,
+ *angle_cosineshift_cost,*angle_cosineshift_theta0;
+ double *angle_cosineshiftexp_umin,*angle_cosineshiftexp_a,
+ *angle_cosineshiftexp_sint,*angle_cosineshiftexp_cost,
+ *angle_cosineshiftexp_theta0;
+ int *angle_cg_cmm_type;
+
+ double *dihedral_charmm_k,*dihedral_charmm_weight;
+ int *dihedral_charmm_multiplicity,*dihedral_charmm_sign;
+ double *dihedral_class2_k1,*dihedral_class2_k2,*dihedral_class2_k3;
+ double *dihedral_class2_phi1,*dihedral_class2_phi2,*dihedral_class2_phi3;
+ double *dihedral_class2_mbt_f1,*dihedral_class2_mbt_f2;
+ double *dihedral_class2_mbt_f3,*dihedral_class2_mbt_r0;
+ double *dihedral_class2_ebt_f1_1,*dihedral_class2_ebt_f2_1;
+ double *dihedral_class2_ebt_f3_1,*dihedral_class2_ebt_r0_1;
+ double *dihedral_class2_ebt_f1_2,*dihedral_class2_ebt_f2_2;
+ double *dihedral_class2_ebt_f3_2,*dihedral_class2_ebt_r0_2;
+ double *dihedral_class2_at_f1_1,*dihedral_class2_at_f2_1;
+ double *dihedral_class2_at_f3_1,*dihedral_class2_at_theta0_1;
+ double *dihedral_class2_at_f1_2,*dihedral_class2_at_f2_2;
+ double *dihedral_class2_at_f3_2,*dihedral_class2_at_theta0_2;
+ double *dihedral_class2_aat_k;
+ double *dihedral_class2_aat_theta0_1,*dihedral_class2_aat_theta0_2;
+ double *dihedral_class2_bb13_k;
+ double *dihedral_class2_bb13_r10,*dihedral_class2_bb13_r30;
+ double *dihedral_harmonic_k;
+ int *dihedral_harmonic_multiplicity,*dihedral_harmonic_sign;
+ double *dihedral_helix_aphi,*dihedral_helix_bphi,*dihedral_helix_cphi;
+ double *dihedral_multi_a1,*dihedral_multi_a2,*dihedral_multi_a3;
+ double *dihedral_multi_a4,*dihedral_multi_a5;
+ double *dihedral_opls_k1,*dihedral_opls_k2;
+ double *dihedral_opls_k3,*dihedral_opls_k4;
+ double *dihedral_cosineshiftexp_umin, *dihedral_cosineshiftexp_a,
+ *dihedral_cosineshiftexp_sint,*dihedral_cosineshiftexp_cost,
+ *dihedral_cosineshiftexp_theta;
+
+ double *improper_class2_k0,*improper_class2_chi0;
+ double *improper_class2_aa_k1,*improper_class2_aa_k2,*improper_class2_aa_k3;
+ double *improper_class2_aa_theta0_1,*improper_class2_aa_theta0_2;
+ double *improper_class2_aa_theta0_3;
+ double *improper_cvff_k;
+ int *improper_cvff_sign,*improper_cvff_multiplicity;
+ double *improper_harmonic_k,*improper_harmonic_chi;
+
+ // atom quantities
+
+ int iatoms,ibonds,iangles,idihedrals,iimpropers;
+
+ double *mass;
+ double *x,*y,*z,*vx,*vy,*vz;
+ double *omegax,*omegay,*omegaz;
+ tagint *tag;
+ int *type,*mask,*image;
+ int *molecule;
+ double *q,*mux,*muy,*muz,*radius,*density,*vfrac,*rmass;
+ double *s0,*x0x,*x0y,*x0z;
+ double *shapex,*shapey,*shapez;
+ double *quatw,*quati,*quatj,*quatk,*angmomx,*angmomy,*angmomz;
+ int *ellipsoid;
+ int *bond_type,*angle_type,*dihedral_type,*improper_type;
+ int *bond_atom1,*bond_atom2;
+ int *angle_atom1,*angle_atom2,*angle_atom3;
+ int *dihedral_atom1,*dihedral_atom2,*dihedral_atom3,*dihedral_atom4;
+ int *improper_atom1,*improper_atom2,*improper_atom3,*improper_atom4;
+
+ // functions
+
+ Data();
+ void stats();
+ void write(FILE *fp, FILE *fp2=NULL);
+
+ void write_atom_angle(FILE *, int, int, int, int);
+ void write_atom_atomic(FILE *, int, int, int, int);
+ void write_atom_bond(FILE *, int, int, int, int);
+ void write_atom_charge(FILE *, int, int, int, int);
+ void write_atom_dipole(FILE *, int, int, int, int);
+ void write_atom_ellipsoid(FILE *, int, int, int, int);
+ void write_atom_full(FILE *, int, int, int, int);
+ void write_atom_molecular(FILE *, int, int, int, int);
+ void write_atom_peri(FILE *, int, int, int, int);
+ void write_atom_sphere(FILE *, int, int, int, int);
+
+ void write_atom_angle_extra(FILE *, int);
+ void write_atom_atomic_extra(FILE *, int);
+ void write_atom_bond_extra(FILE *, int);
+ void write_atom_charge_extra(FILE *, int);
+ void write_atom_dipole_extra(FILE *, int);
+ void write_atom_ellipsoid_extra(FILE *, int);
+ void write_atom_full_extra(FILE *, int);
+ void write_atom_molecular_extra(FILE *, int);
+ void write_atom_peri_extra(FILE *, int);
+ void write_atom_sphere_extra(FILE *, int);
+
+ void write_vel_angle(FILE *, int);
+ void write_vel_atomic(FILE *, int);
+ void write_vel_bond(FILE *, int);
+ void write_vel_charge(FILE *, int);
+ void write_vel_dipole(FILE *, int);
+ void write_vel_ellipsoid(FILE *, int);
+ void write_vel_full(FILE *, int);
+ void write_vel_molecular(FILE *, int);
+ void write_vel_peri(FILE *, int);
+ void write_vel_sphere(FILE *, int);
+
+ void write_vel_angle_extra(FILE *, int);
+ void write_vel_atomic_extra(FILE *, int);
+ void write_vel_bond_extra(FILE *, int);
+ void write_vel_charge_extra(FILE *, int);
+ void write_vel_dipole_extra(FILE *, int);
+ void write_vel_ellipsoid_extra(FILE *, int);
+ void write_vel_full_extra(FILE *, int);
+ void write_vel_molecular_extra(FILE *, int);
+ void write_vel_peri_extra(FILE *, int);
+ void write_vel_sphere_extra(FILE *, int);
+};
+
+// ---------------------------------------------------------------------
+// function prototypes
+// ---------------------------------------------------------------------
+
+void header(FILE *, Data &);
+void set_style(char *, Data &, int);
+void groups(FILE *);
+void type_arrays(FILE *, Data &);
+void force_fields(FILE *, Data &);
+void modify(FILE *);
+void pair(FILE *fp, Data &data, char *style, int flag);
+void bond(FILE *fp, Data &data);
+void angle(FILE *fp, Data &data);
+void dihedral(FILE *fp, Data &data);
+void improper(FILE *fp, Data &data);
+int atom(double *, Data &data);
+
+void allocate_angle(Data &data);
+void allocate_atomic(Data &data);
+void allocate_bond(Data &data);
+void allocate_charge(Data &data);
+void allocate_dipole(Data &data);
+void allocate_ellipsoid(Data &data);
+void allocate_full(Data &data);
+void allocate_molecular(Data &data);
+void allocate_peri(Data &data);
+void allocate_sphere(Data &data);
+
+int atom_angle(double *, Data &, int);
+int atom_atomic(double *, Data &, int);
+int atom_bond(double *, Data &, int);
+int atom_charge(double *, Data &, int);
+int atom_dipole(double *, Data &, int);
+int atom_ellipsoid(double *, Data &, int);
+int atom_full(double *, Data &, int);
+int atom_molecular(double *, Data &, int);
+int atom_peri(double *, Data &, int);
+int atom_sphere(double *, Data &, int);
+
+void strip_suffix(char *);
+
+int read_int(FILE *fp);
+double read_double(FILE *fp);
+char *read_char(FILE *fp);
+bigint read_bigint(FILE *fp);
+
+// ---------------------------------------------------------------------
+// main program
+// ---------------------------------------------------------------------
+
+int main (int narg, char **arg)
+{
+ // process command-line args
+
+ int iarg = 1;
+ if (strcmp(arg[iarg],"-h") == 0) {
+ printf("Syntax: restart2data switch arg ... "
+ "restart-file data-file (input-file)\n");
+ printf(" restart-file and data-file are mandatory");
+ printf(" input-file is optional");
+ printf(" if specified it will contain LAMMPS input script commands");
+ printf(" for mass and force field info");
+ printf(" only a few force field styles support this option");
+ return 0;
+ }
+
+ if ((narg-iarg != 2) && (narg-iarg != 3)) {
+ printf("Syntax: restart2data switch arg ... "
+ "restart-file data-file (input-file)\n");
+ return 1;
+ }
+
+ char *restartfile = arg[iarg];
+ char *datafile = arg[iarg+1];
+ char *inputfile = NULL;
+ if (narg-iarg == 3) inputfile = arg[iarg+2];
+
+ // if restart file contains '%', file = filename with % replaced by "base"
+ // else file = single file
+ // open single restart file or base file for multiproc case
+
+ printf("Reading restart file ...\n");
+
+ char *ptr;
+ FILE *fp;
+
+ int multiproc = 0;
+ if ( (ptr = strchr(restartfile,'%')) ) {
+ multiproc = 1;
+ char *basefile = new char[strlen(restartfile) + 16];
+ *ptr = '\0';
+ sprintf(basefile,"%s%s%s",restartfile,"base",ptr+1);
+ fp = fopen(basefile,"rb");
+ if (fp == NULL) {
+ printf("ERROR: Cannot open restart file %s\n",basefile);
+ return 1;
+ }
+ } else {
+ fp = fopen(restartfile,"rb");
+ if (fp == NULL) {
+ printf("ERROR: Cannot open restart file %s\n",restartfile);
+ return 1;
+ }
+ }
+
+ // read beginning of restart file
+
+ Data data;
+
+ header(fp,data);
+ if (data.size_smallint != sizeof(int) ||
+ data.size_tagint != sizeof(tagint) ||
+ data.size_bigint != sizeof(bigint)) {
+ printf("ERROR: Data type sizes in restart file "
+ "are incompatible with restart2data.cpp\n");
+ return 1;
+ }
+
+ groups(fp);
+ type_arrays(fp,data);
+ force_fields(fp,data);
+ modify(fp);
+
+ // read atoms from single or multiple restart files
+
+ double *buf = NULL;
+ int n,m;
+ int maxbuf = 0;
+ data.iatoms = data.ibonds = data.iangles =
+ data.idihedrals = data.iimpropers = 0;
+
+ for (int iproc = 0; iproc < data.nprocs; iproc++) {
+ if (multiproc) {
+ fclose(fp);
+ char *procfile;
+ sprintf(procfile,"%s%d%s",restartfile,iproc,ptr+1);
+ fp = fopen(procfile,"rb");
+ if (fp == NULL) {
+ printf("ERROR: Cannot open restart file %s\n",procfile);
+ return 1;
+ }
+ }
+ n = read_int(fp);
+
+ if (n > maxbuf) {
+ maxbuf = n;
+ delete [] buf;
+ buf = new double[maxbuf];
+ }
+
+ fread(buf,sizeof(double),n,fp);
+
+ m = 0;
+ while (m < n) m += atom(&buf[m],data);
+ }
+
+ fclose(fp);
+
+ // print out stats
+
+ data.stats();
+
+ // write out data file and no input file
+
+ if (!inputfile) {
+ printf("Writing data file ...\n");
+ fp = fopen(datafile,"w");
+ if (fp == NULL) {
+ printf("ERROR: Cannot open data file %s\n",datafile);
+ return 1;
+ }
+ data.write(fp);
+ fclose(fp);
+
+ // write out data file and input file
+
+ } else {
+ printf("Writing data file ...\n");
+ fp = fopen(datafile,"w");
+ if (fp == NULL) {
+ printf("ERROR: Cannot open data file %s\n",datafile);
+ return 1;
+ }
+ printf("Writing input file ...\n");
+ FILE *fp2 = fopen(inputfile,"w");
+ if (fp2 == NULL) {
+ printf("ERROR: Cannot open input file %s\n",inputfile);
+ return 1;
+ }
+
+ data.write(fp,fp2);
+ fclose(fp);
+ fclose(fp2);
+ }
+
+ return 0;
+}
+
+// ---------------------------------------------------------------------
+// read header of restart file
+// ---------------------------------------------------------------------
+
+void header(FILE *fp, Data &data)
+{
+ char *version = "19 Aug 2011";
+
+ data.triclinic = 0;
+
+ int flag;
+ flag = read_int(fp);
+
+ while (flag >= 0) {
+
+ if (flag == VERSION) {
+ data.version = read_char(fp);
+ if (strcmp(version,data.version) != 0) {
+ char *str = "Restart file version does not match restart2data version";
+ printf("WARNING %s\n",str);
+ printf(" restart2data version = %s\n",version);
+ }
+ }
+ else if (flag == SMALLINT) data.size_smallint = read_int(fp);
+ else if (flag == TAGINT) data.size_tagint = read_int(fp);
+ else if (flag == BIGINT) data.size_bigint = read_int(fp);
+ else if (flag == UNITS) data.unit_style = read_char(fp);
+ else if (flag == NTIMESTEP) data.ntimestep = read_bigint(fp);
+ else if (flag == DIMENSION) data.dimension = read_int(fp);
+ else if (flag == NPROCS) data.nprocs = read_int(fp);
+ else if (flag == PROCGRID_0) data.px = read_int(fp);
+ else if (flag == PROCGRID_1) data.py = read_int(fp);
+ else if (flag == PROCGRID_2) data.pz = read_int(fp);
+ else if (flag == NEWTON_PAIR) data.newton_pair = read_int(fp);
+ else if (flag == NEWTON_BOND) data.newton_bond = read_int(fp);
+ else if (flag == XPERIODIC) data.xperiodic = read_int(fp);
+ else if (flag == YPERIODIC) data.yperiodic = read_int(fp);
+ else if (flag == ZPERIODIC) data.zperiodic = read_int(fp);
+ else if (flag == BOUNDARY_00) data.boundary[0][0] = read_int(fp);
+ else if (flag == BOUNDARY_01) data.boundary[0][1] = read_int(fp);
+ else if (flag == BOUNDARY_10) data.boundary[1][0] = read_int(fp);
+ else if (flag == BOUNDARY_11) data.boundary[1][1] = read_int(fp);
+ else if (flag == BOUNDARY_20) data.boundary[2][0] = read_int(fp);
+ else if (flag == BOUNDARY_21) data.boundary[2][1] = read_int(fp);
+
+ // if atom_style = hybrid:
+ // set data_style_hybrid to # of sub-styles
+ // read additional sub-class arguments
+ // set sub-styles to 1 to N
+
+ else if (flag == ATOM_STYLE) {
+ data.style_angle = data.style_atomic = data.style_bond =
+ data.style_charge = data.style_dipole =
+ data.style_ellipsoid = data.style_full = data.style_hybrid =
+ data.style_molecular = data.style_peri = data.style_sphere = 0;
+
+ data.atom_style = read_char(fp);
+ strip_suffix(data.atom_style);
+ set_style(data.atom_style,data,1);
+
+ if (strcmp(data.atom_style,"hybrid") == 0) {
+ int nwords = read_int(fp);
+ set_style(data.atom_style,data,nwords);
+ char *substyle;
+ for (int i = 1; i <= nwords; i++) {
+ substyle = read_char(fp);
+ set_style(substyle,data,i);
+ }
+ }
+ }
+
+ else if (flag == NATOMS) data.natoms = read_bigint(fp);
+ else if (flag == NTYPES) data.ntypes = read_int(fp);
+ else if (flag == NBONDS) data.nbonds = read_bigint(fp);
+ else if (flag == NBONDTYPES) data.nbondtypes = read_int(fp);
+ else if (flag == BOND_PER_ATOM) data.bond_per_atom = read_int(fp);
+ else if (flag == NANGLES) data.nangles = read_bigint(fp);
+ else if (flag == NANGLETYPES) data.nangletypes = read_int(fp);
+ else if (flag == ANGLE_PER_ATOM) data.angle_per_atom = read_int(fp);
+ else if (flag == NDIHEDRALS) data.ndihedrals = read_bigint(fp);
+ else if (flag == NDIHEDRALTYPES) data.ndihedraltypes = read_int(fp);
+ else if (flag == DIHEDRAL_PER_ATOM) data.dihedral_per_atom = read_int(fp);
+ else if (flag == NIMPROPERS) data.nimpropers = read_bigint(fp);
+ else if (flag == NIMPROPERTYPES) data.nimpropertypes = read_int(fp);
+ else if (flag == IMPROPER_PER_ATOM) data.improper_per_atom = read_int(fp);
+ else if (flag == BOXLO_0) data.xlo = read_double(fp);
+ else if (flag == BOXHI_0) data.xhi = read_double(fp);
+ else if (flag == BOXLO_1) data.ylo = read_double(fp);
+ else if (flag == BOXHI_1) data.yhi = read_double(fp);
+ else if (flag == BOXLO_2) data.zlo = read_double(fp);
+ else if (flag == BOXHI_2) data.zhi = read_double(fp);
+ else if (flag == SPECIAL_LJ_1) data.special_lj[1] = read_double(fp);
+ else if (flag == SPECIAL_LJ_2) data.special_lj[2] = read_double(fp);
+ else if (flag == SPECIAL_LJ_3) data.special_lj[3] = read_double(fp);
+ else if (flag == SPECIAL_COUL_1) data.special_coul[1] = read_double(fp);
+ else if (flag == SPECIAL_COUL_2) data.special_coul[2] = read_double(fp);
+ else if (flag == SPECIAL_COUL_3) data.special_coul[3] = read_double(fp);
+ else if (flag == XY) {
+ data.triclinic = 1;
+ data.xy = read_double(fp);
+ } else if (flag == XZ) {
+ data.triclinic = 1;
+ data.xz = read_double(fp);
+ } else if (flag == YZ) {
+ data.triclinic = 1;
+ data.yz = read_double(fp);
+ } else {
+ printf("ERROR: Invalid flag in header section of restart file %d\n",
+ flag);
+ exit(1);
+ }
+
+ flag = read_int(fp);
+ }
+}
+
+// ---------------------------------------------------------------------
+// set atom style to flag
+// ---------------------------------------------------------------------
+
+void set_style(char *name, Data &data, int flag)
+{
+ if (strcmp(name,"angle") == 0) data.style_angle = flag;
+ else if (strcmp(name,"atomic") == 0) data.style_atomic = flag;
+ else if (strcmp(name,"bond") == 0) data.style_bond = flag;
+ else if (strcmp(name,"charge") == 0) data.style_charge = flag;
+ else if (strcmp(name,"dipole") == 0) data.style_dipole = flag;
+ else if (strcmp(name,"ellipsoid") == 0) data.style_ellipsoid = flag;
+ else if (strcmp(name,"full") == 0) data.style_full = flag;
+ else if (strcmp(name,"hybrid") == 0) data.style_hybrid = flag;
+ else if (strcmp(name,"molecular") == 0) data.style_molecular = flag;
+ else if (strcmp(name,"peri") == 0) data.style_peri = flag;
+ else if (strcmp(name,"sphere") == 0) data.style_sphere = flag;
+ else {
+ printf("ERROR: Unknown atom style %s\n",name);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// read group info from restart file, just ignore it
+// ---------------------------------------------------------------------
+
+void groups(FILE *fp)
+{
+ int ngroup = read_int(fp);
+
+ char *name;
+
+ // use count to not change restart format with deleted groups
+ // remove this on next major release
+
+ int count = 0;
+ for (int i = 0; i < MAX_GROUP; i++) {
+ name = read_char(fp);
+ delete [] name;
+ count++;
+ if (count == ngroup) break;
+ }
+}
+
+// ---------------------------------------------------------------------
+// read type arrays from restart file
+// ---------------------------------------------------------------------
+
+void type_arrays(FILE *fp, Data &data)
+{
+ data.mass = NULL;
+
+ int flag;
+ flag = read_int(fp);
+
+ while (flag >= 0) {
+
+ if (flag == MASS) {
+ data.mass = new double[data.ntypes+1];
+ fread(&data.mass[1],sizeof(double),data.ntypes,fp);
+ } else {
+ printf("ERROR: Invalid flag in type arrays section of restart file %d\n",
+ flag);
+ exit(1);
+ }
+
+ flag = read_int(fp);
+ }
+}
+
+// ---------------------------------------------------------------------
+// read force-field info from restart file
+// ---------------------------------------------------------------------
+
+void force_fields(FILE *fp, Data &data)
+{
+ data.pair_style = data.bond_style = data.angle_style =
+ data.dihedral_style = data.improper_style = NULL;
+
+ int flag;
+ flag = read_int(fp);
+
+ while (flag >= 0) {
+
+ if (flag == PAIR) {
+ data.pair_style = read_char(fp);
+ strip_suffix(data.pair_style);
+ pair(fp,data,data.pair_style,1);
+ } else if (flag == BOND) {
+ data.bond_style = read_char(fp);
+ strip_suffix(data.bond_style);
+ bond(fp,data);
+ } else if (flag == ANGLE) {
+ data.angle_style = read_char(fp);
+ strip_suffix(data.angle_style);
+ angle(fp,data);
+ } else if (flag == DIHEDRAL) {
+ data.dihedral_style = read_char(fp);
+ strip_suffix(data.dihedral_style);
+ dihedral(fp,data);
+ } else if (flag == IMPROPER) {
+ data.improper_style = read_char(fp);
+ strip_suffix(data.improper_style);
+ improper(fp,data);
+ } else {
+ printf("ERROR: Invalid flag in force fields section of restart file %d\n",
+ flag);
+ exit(1);
+ }
+
+ flag = read_int(fp);
+ }
+}
+
+// ---------------------------------------------------------------------
+// read fix info from restart file, just ignore it
+// ---------------------------------------------------------------------
+
+void modify(FILE *fp)
+{
+ char *buf;
+ int n;
+
+ // nfix = # of fix entries with state
+
+ int nfix = read_int(fp);
+
+ // read each entry with id string, style string, chunk of data
+
+ for (int i = 0; i < nfix; i++) {
+ buf = read_char(fp); delete [] buf;
+ buf = read_char(fp); delete [] buf;
+ buf = read_char(fp); delete [] buf;
+ }
+
+ // nfix = # of fix entries with peratom info
+
+ int nfix_peratom = read_int(fp);
+
+ // read each entry with id string, style string, maxsize of one atom data
+
+ for (int i = 0; i < nfix_peratom; i++) {
+ buf = read_char(fp); delete [] buf;
+ buf = read_char(fp); delete [] buf;
+ n = read_int(fp);
+ }
+}
+
+// ---------------------------------------------------------------------
+// read atom info from restart file and store in data struct
+// ---------------------------------------------------------------------
+
+int atom(double *buf, Data &data)
+{
+ // allocate per-atom arrays
+
+ if (data.iatoms == 0) {
+
+ // common to all atom styles
+
+ data.x = new double[data.natoms];
+ data.y = new double[data.natoms];
+ data.z = new double[data.natoms];
+ data.tag = new int[data.natoms];
+ data.type = new int[data.natoms];
+ data.mask = new int[data.natoms];
+ data.image = new int[data.natoms];
+ data.vx = new double[data.natoms];
+ data.vy = new double[data.natoms];
+ data.vz = new double[data.natoms];
+
+ // style-specific arrays
+ // don't worry about re-allocating if style = hybrid
+
+ if (data.style_angle) allocate_angle(data);
+ if (data.style_atomic) allocate_atomic(data);
+ if (data.style_bond) allocate_bond(data);
+ if (data.style_charge) allocate_charge(data);
+ if (data.style_dipole) allocate_dipole(data);
+ if (data.style_ellipsoid) allocate_ellipsoid(data);
+ if (data.style_full) allocate_full(data);
+ if (data.style_molecular) allocate_molecular(data);
+ if (data.style_peri) allocate_peri(data);
+ if (data.style_sphere) allocate_sphere(data);
+ }
+
+ // read atom quantities from buf
+ // if hybrid, loop over all sub-styles in order listed
+ // if hybrid, loop index k will match style setting to insure correct order
+
+ int nloop = 1;
+ if (data.style_hybrid) nloop = data.style_hybrid;
+
+ int iatoms = data.iatoms;
+ int m = 0;
+ for (int k = 1; k <= nloop; k++) {
+ if (k == data.style_angle) m += atom_angle(&buf[m],data,iatoms);
+ if (k == data.style_atomic) m += atom_atomic(&buf[m],data,iatoms);
+ if (k == data.style_bond) m += atom_bond(&buf[m],data,iatoms);
+ if (k == data.style_charge) m += atom_charge(&buf[m],data,iatoms);
+ if (k == data.style_dipole) m += atom_dipole(&buf[m],data,iatoms);
+ if (k == data.style_ellipsoid) m += atom_ellipsoid(&buf[m],data,iatoms);
+ if (k == data.style_full) m += atom_full(&buf[m],data,iatoms);
+ if (k == data.style_molecular) m += atom_molecular(&buf[m],data,iatoms);
+ if (k == data.style_peri) m += atom_peri(&buf[m],data,iatoms);
+ if (k == data.style_sphere) m += atom_sphere(&buf[m],data,iatoms);
+ }
+
+ data.iatoms++;
+ m = static_cast<int> (buf[0]);
+ return m;
+}
+
+// ---------------------------------------------------------------------
+// read one atom's info from buffer
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+int atom_angle(double *buf, Data &data, int iatoms)
+{
+ int type,atom1,atom2,atom3;
+
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+ int n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] < atom1 ) {
+ data.bond_type[data.ibonds] = type;
+ data.bond_atom1[data.ibonds] = data.tag[iatoms];
+ data.bond_atom2[data.ibonds] = atom1;
+ data.ibonds++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.angle_type[data.iangles] = type;
+ data.angle_atom1[data.iangles] = atom1;
+ data.angle_atom2[data.iangles] = atom2;
+ data.angle_atom3[data.iangles] = atom3;
+ data.iangles++;
+ }
+ }
+
+ return m;
+}
+
+int atom_atomic(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ return m;
+}
+
+int atom_bond(double *buf, Data &data, int iatoms)
+{
+ int type,atom1;
+
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+ int n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] < atom1) {
+ data.bond_type[data.ibonds] = type;
+ data.bond_atom1[data.ibonds] = data.tag[iatoms];
+ data.bond_atom2[data.ibonds] = atom1;
+ data.ibonds++;
+ }
+ }
+
+ return m;
+}
+
+int atom_charge(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.q[iatoms] = buf[m++];
+
+ return m;
+}
+
+int atom_dipole(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.q[iatoms] = buf[m++];
+ data.mux[iatoms] = buf[m++];
+ data.muy[iatoms] = buf[m++];
+ data.muz[iatoms] = buf[m++];
+
+ return m;
+}
+
+int atom_ellipsoid(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.rmass[iatoms] = buf[m++];
+ data.angmomx[iatoms] = buf[m++];
+ data.angmomy[iatoms] = buf[m++];
+ data.angmomz[iatoms] = buf[m++];
+ data.ellipsoid[iatoms] = static_cast<int> (buf[m++]);
+
+ if (data.ellipsoid[iatoms]) {
+ data.nellipsoids++;
+ data.shapex[iatoms] = buf[m++];
+ data.shapey[iatoms] = buf[m++];
+ data.shapez[iatoms] = buf[m++];
+ data.quatw[iatoms] = buf[m++];
+ data.quati[iatoms] = buf[m++];
+ data.quatj[iatoms] = buf[m++];
+ data.quatk[iatoms] = buf[m++];
+ data.density[iatoms] = data.rmass[iatoms] /
+ (4.0*PI/3.0 *
+ data.shapex[iatoms]*data.shapey[iatoms]*data.shapez[iatoms]);
+ } else data.density[iatoms] = data.rmass[iatoms];
+
+ return m;
+}
+
+int atom_sphere(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.radius[iatoms] = buf[m++];
+ data.rmass[iatoms] = buf[m++];
+ if (data.radius[iatoms] == 0.0) data.density[iatoms] = data.rmass[iatoms];
+ else
+ data.density[iatoms] = data.rmass[iatoms] /
+ (4.0*PI/3.0 *
+ data.radius[iatoms]*data.radius[iatoms]*data.radius[iatoms]);
+ data.omegax[iatoms] = buf[m++];
+ data.omegay[iatoms] = buf[m++];
+ data.omegaz[iatoms] = buf[m++];
+
+ return m;
+}
+
+int atom_full(double *buf, Data &data, int iatoms)
+{
+ int type,atom1,atom2,atom3,atom4;
+
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.q[iatoms] = buf[m++];
+ data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+ int n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] < atom1) {
+ data.bond_type[data.ibonds] = type;
+ data.bond_atom1[data.ibonds] = data.tag[iatoms];
+ data.bond_atom2[data.ibonds] = atom1;
+ data.ibonds++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.angle_type[data.iangles] = type;
+ data.angle_atom1[data.iangles] = atom1;
+ data.angle_atom2[data.iangles] = atom2;
+ data.angle_atom3[data.iangles] = atom3;
+ data.iangles++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ atom4 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.dihedral_type[data.idihedrals] = type;
+ data.dihedral_atom1[data.idihedrals] = atom1;
+ data.dihedral_atom2[data.idihedrals] = atom2;
+ data.dihedral_atom3[data.idihedrals] = atom3;
+ data.dihedral_atom4[data.idihedrals] = atom4;
+ data.idihedrals++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ atom4 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.improper_type[data.iimpropers] = type;
+ data.improper_atom1[data.iimpropers] = atom1;
+ data.improper_atom2[data.iimpropers] = atom2;
+ data.improper_atom3[data.iimpropers] = atom3;
+ data.improper_atom4[data.iimpropers] = atom4;
+ data.iimpropers++;
+ }
+ }
+
+ return m;
+}
+
+int atom_molecular(double *buf, Data &data, int iatoms)
+{
+ int type,atom1,atom2,atom3,atom4;
+
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+ int n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] < atom1) {
+ data.bond_type[data.ibonds] = type;
+ data.bond_atom1[data.ibonds] = data.tag[iatoms];
+ data.bond_atom2[data.ibonds] = atom1;
+ data.ibonds++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.angle_type[data.iangles] = type;
+ data.angle_atom1[data.iangles] = atom1;
+ data.angle_atom2[data.iangles] = atom2;
+ data.angle_atom3[data.iangles] = atom3;
+ data.iangles++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ atom4 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.dihedral_type[data.idihedrals] = type;
+ data.dihedral_atom1[data.idihedrals] = atom1;
+ data.dihedral_atom2[data.idihedrals] = atom2;
+ data.dihedral_atom3[data.idihedrals] = atom3;
+ data.dihedral_atom4[data.idihedrals] = atom4;
+ data.idihedrals++;
+ }
+ }
+
+ n = static_cast<int> (buf[m++]);
+ for (int k = 0; k < n; k++) {
+ type = static_cast<int> (buf[m++]);
+ atom1 = static_cast<int> (buf[m++]);
+ atom2 = static_cast<int> (buf[m++]);
+ atom3 = static_cast<int> (buf[m++]);
+ atom4 = static_cast<int> (buf[m++]);
+ if (data.newton_bond || data.tag[iatoms] == atom2) {
+ data.improper_type[data.iimpropers] = type;
+ data.improper_atom1[data.iimpropers] = atom1;
+ data.improper_atom2[data.iimpropers] = atom2;
+ data.improper_atom3[data.iimpropers] = atom3;
+ data.improper_atom4[data.iimpropers] = atom4;
+ data.iimpropers++;
+ }
+ }
+
+ return m;
+}
+
+int atom_peri(double *buf, Data &data, int iatoms)
+{
+ int m = 1;
+ data.x[iatoms] = buf[m++];
+ data.y[iatoms] = buf[m++];
+ data.z[iatoms] = buf[m++];
+ data.tag[iatoms] = static_cast<int> (buf[m++]);
+ data.type[iatoms] = static_cast<int> (buf[m++]);
+ data.mask[iatoms] = static_cast<int> (buf[m++]);
+ data.image[iatoms] = static_cast<int> (buf[m++]);
+ data.vx[iatoms] = buf[m++];
+ data.vy[iatoms] = buf[m++];
+ data.vz[iatoms] = buf[m++];
+
+ data.vfrac[iatoms] = buf[m++];
+ data.rmass[iatoms] = buf[m++];
+ data.s0[iatoms] = buf[m++];
+ data.x0x[iatoms] = buf[m++];
+ data.x0y[iatoms] = buf[m++];
+ data.x0z[iatoms] = buf[m++];
+
+ return m;
+}
+
+// ---------------------------------------------------------------------
+// per-atom memory allocation routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void allocate_angle(Data &data)
+{
+ data.molecule = new int[data.natoms];
+ data.bond_type = new int[data.nbonds];
+ data.bond_atom1 = new int[data.nbonds];
+ data.bond_atom2 = new int[data.nbonds];
+ data.angle_type = new int[data.nangles];
+ data.angle_atom1 = new int[data.nangles];
+ data.angle_atom2 = new int[data.nangles];
+ data.angle_atom3 = new int[data.nangles];
+}
+
+void allocate_atomic(Data &data) {}
+
+void allocate_bond(Data &data)
+{
+ data.molecule = new int[data.natoms];
+ data.bond_type = new int[data.nbonds];
+ data.bond_atom1 = new int[data.nbonds];
+ data.bond_atom2 = new int[data.nbonds];
+}
+
+void allocate_charge(Data &data)
+{
+ data.q = new double[data.natoms];
+}
+
+void allocate_dipole(Data &data)
+{
+ data.q = new double[data.natoms];
+ data.mux = new double[data.natoms];
+ data.muy = new double[data.natoms];
+ data.muz = new double[data.natoms];
+}
+
+void allocate_full(Data &data)
+{
+ data.q = new double[data.natoms];
+ data.molecule = new int[data.natoms];
+ data.bond_type = new int[data.nbonds];
+ data.bond_atom1 = new int[data.nbonds];
+ data.bond_atom2 = new int[data.nbonds];
+ data.angle_type = new int[data.nangles];
+ data.angle_atom1 = new int[data.nangles];
+ data.angle_atom2 = new int[data.nangles];
+ data.angle_atom3 = new int[data.nangles];
+ data.dihedral_type = new int[data.ndihedrals];
+ data.dihedral_atom1 = new int[data.ndihedrals];
+ data.dihedral_atom2 = new int[data.ndihedrals];
+ data.dihedral_atom3 = new int[data.ndihedrals];
+ data.dihedral_atom4 = new int[data.ndihedrals];
+ data.improper_type = new int[data.nimpropers];
+ data.improper_atom1 = new int[data.nimpropers];
+ data.improper_atom2 = new int[data.nimpropers];
+ data.improper_atom3 = new int[data.nimpropers];
+ data.improper_atom4 = new int[data.nimpropers];
+}
+
+void allocate_ellipsoid(Data &data)
+{
+ data.rmass = new double[data.natoms];
+ data.density = new double[data.natoms];
+ data.angmomx = new double[data.natoms];
+ data.angmomy = new double[data.natoms];
+ data.angmomz = new double[data.natoms];
+ data.ellipsoid = new int[data.natoms];
+ data.quatw = new double[data.natoms];
+ data.shapex = new double[data.natoms];
+ data.shapey = new double[data.natoms];
+ data.shapez = new double[data.natoms];
+ data.quati = new double[data.natoms];
+ data.quatj = new double[data.natoms];
+ data.quatk = new double[data.natoms];
+}
+
+void allocate_sphere(Data &data)
+{
+ data.radius = new double[data.natoms];
+ data.rmass = new double[data.natoms];
+ data.density = new double[data.natoms];
+ data.omegax = new double[data.natoms];
+ data.omegay = new double[data.natoms];
+ data.omegaz = new double[data.natoms];
+}
+
+void allocate_molecular(Data &data)
+{
+ data.molecule = new int[data.natoms];
+ data.bond_type = new int[data.nbonds];
+ data.bond_atom1 = new int[data.nbonds];
+ data.bond_atom2 = new int[data.nbonds];
+ data.angle_type = new int[data.nangles];
+ data.angle_atom1 = new int[data.nangles];
+ data.angle_atom2 = new int[data.nangles];
+ data.angle_atom3 = new int[data.nangles];
+ data.dihedral_type = new int[data.ndihedrals];
+ data.dihedral_atom1 = new int[data.ndihedrals];
+ data.dihedral_atom2 = new int[data.ndihedrals];
+ data.dihedral_atom3 = new int[data.ndihedrals];
+ data.dihedral_atom4 = new int[data.ndihedrals];
+ data.improper_type = new int[data.nimpropers];
+ data.improper_atom1 = new int[data.nimpropers];
+ data.improper_atom2 = new int[data.nimpropers];
+ data.improper_atom3 = new int[data.nimpropers];
+ data.improper_atom4 = new int[data.nimpropers];
+}
+
+void allocate_peri(Data &data)
+{
+ data.vfrac = new double[data.natoms];
+ data.rmass = new double[data.natoms];
+ data.s0 = new double[data.natoms];
+ data.x0x = new double[data.natoms];
+ data.x0y = new double[data.natoms];
+ data.x0z = new double[data.natoms];
+}
+
+// ---------------------------------------------------------------------
+// pair coeffs
+// one section for each pair style
+// flag = 1, read all coeff info and allocation arrays
+// flag = 0, just read global settings (when called recursively by hybrid)
+// ---------------------------------------------------------------------
+
+void pair(FILE *fp, Data &data, char *style, int flag)
+{
+ int i,j,m;
+ int itmp;
+ double rtmp;
+
+ if (strcmp(style,"none") == 0) {
+
+ } else if (strcmp(style,"adp") == 0) {
+ } else if (strcmp(style,"airebo") == 0) {
+
+ } else if (strcmp(style,"born/coul/long") == 0) {
+ double cut_lj_global = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_born_A = new double[data.ntypes+1];
+ data.pair_born_rho = new double[data.ntypes+1];
+ data.pair_born_sigma = new double[data.ntypes+1];
+ data.pair_born_C = new double[data.ntypes+1];
+ data.pair_born_D = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_born_A[i] = read_double(fp);
+ data.pair_born_rho[i] = read_double(fp);
+ data.pair_born_sigma[i] = read_double(fp);
+ data.pair_born_C[i] = read_double(fp);
+ data.pair_born_D[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ } else {
+ double born_A = read_double(fp);
+ double born_rho = read_double(fp);
+ double born_sigma = read_double(fp);
+ double born_C = read_double(fp);
+ double born_D = read_double(fp);
+ double cut_lj = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"buck") == 0) ||
+ (strcmp(style,"buck/coul/cut") == 0) ||
+ (strcmp(style,"buck/coul/long") == 0) ||
+ (strcmp(style,"buck/coul") == 0)) {
+
+ if (strcmp(style,"buck") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"buck/coul/cut") == 0) {
+ m = 1;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"buck/coul/long") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"buck/coul") == 0) {
+ m = 0;
+ double cut_buck_global = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ int ewald_order = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ data.pair_buck_A = new double[data.ntypes+1];
+ data.pair_buck_rho = new double[data.ntypes+1];
+ data.pair_buck_C = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_buck_A[i] = read_double(fp);
+ data.pair_buck_rho[i] = read_double(fp);
+ data.pair_buck_C[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ } else {
+ double buck_A = read_double(fp);
+ double buck_rho = read_double(fp);
+ double buck_C = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"colloid") == 0) {
+
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_colloid_A12 = new double[data.ntypes+1];
+ data.pair_colloid_sigma = new double[data.ntypes+1];
+ data.pair_colloid_d1 = new double[data.ntypes+1];
+ data.pair_colloid_d2 = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_colloid_A12[i] = read_double(fp);
+ data.pair_colloid_sigma[i] = read_double(fp);
+ data.pair_colloid_d1[i] = read_double(fp);
+ data.pair_colloid_d2[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ } else {
+ double colloid_A12 = read_double(fp);
+ double colloid_sigma = read_double(fp);
+ double colloid_d1 = read_double(fp);
+ double colloid_d2 = read_double(fp);
+ double cut_lj = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"coul/cut") == 0) ||
+ (strcmp(style,"coul/debye") == 0) ||
+ (strcmp(style,"coul/long") == 0)) {
+
+ if (strcmp(style,"coul/cut") == 0) {
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"coul/debye") == 0) {
+ m = 1;
+ double cut_coul = read_double(fp);
+ double kappa = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"coul/long") == 0) {
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ double cut_coul = read_double(fp);
+ } else {
+ double cut_coul = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"dipole/cut") == 0) {
+
+ double cut_lj_global = read_double(fp);
+ double cut_coul_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_dipole_epsilon = new double[data.ntypes+1];
+ data.pair_dipole_sigma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_dipole_epsilon[i] = read_double(fp);
+ data.pair_dipole_sigma[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul = read_double(fp);
+ } else {
+ double dipole_epsilon = read_double(fp);
+ double dipole_sigma = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"dpd") == 0) {
+
+ double temperature = read_double(fp);
+ double cut_global = read_double(fp);
+ int seed = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_dpd_a0 = new double[data.ntypes+1];
+ data.pair_dpd_gamma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_dpd_a0[i] = read_double(fp);
+ data.pair_dpd_gamma[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double dpd_a0 = read_double(fp);
+ double dpd_gamma = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"dpd/tstat") == 0) {
+
+ double tstart = read_double(fp);
+ double tstop = read_double(fp);
+ double cut_global = read_double(fp);
+ int seed = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_dpd_a0 = new double[data.ntypes+1];
+ data.pair_dpd_gamma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_dpd_a0[i] = read_double(fp);
+ data.pair_dpd_gamma[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double dpd_a0 = read_double(fp);
+ double dpd_gamma = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"eam") == 0) {
+ } else if (strcmp(style,"eam/alloy") == 0) {
+ } else if (strcmp(style,"eam/fs") == 0) {
+ } else if (strcmp(style,"eim") == 0) {
+
+ } else if (strcmp(style,"eff/cut") == 0) {
+
+ double cut_coul = read_double(fp);
+ int limit_size_flag = read_int(fp);
+ int flexible_pressure_flag = read_int(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ double cut = read_double(fp);
+ } else {
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"gayberne") == 0) {
+
+ double gamma = read_double(fp);
+ double upsilon = read_double(fp);
+ double mu = read_double(fp);
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_gb_epsilon = new double[data.ntypes+1];
+ data.pair_gb_sigma = new double[data.ntypes+1];
+ data.pair_gb_epsa = new double[data.ntypes+1];
+ data.pair_gb_epsb = new double[data.ntypes+1];
+ data.pair_gb_epsc = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++) {
+ itmp = read_int(fp);
+ if (itmp) {
+ data.pair_gb_epsa[i] = read_double(fp);
+ data.pair_gb_epsb[i] = read_double(fp);
+ data.pair_gb_epsc[i] = read_double(fp);
+ data.pair_gb_epsa[i] = pow(data.pair_gb_epsa[i],-mu);
+ data.pair_gb_epsb[i] = pow(data.pair_gb_epsb[i],-mu);
+ data.pair_gb_epsc[i] = pow(data.pair_gb_epsc[i],-mu);
+ }
+
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_gb_epsilon[i] = read_double(fp);
+ data.pair_gb_sigma[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double gb_epsilon = read_double(fp);
+ double gb_sigma = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"gran/hooke") == 0) ||
+ (strcmp(style,"gran/hooke/history") == 0) ||
+ (strcmp(style,"gran/hertz/history") == 0)) {
+
+ double kn = read_double(fp);
+ double kt = read_double(fp);
+ double gamman = read_double(fp);
+ double gammat = read_double(fp);
+ double xmu = read_double(fp);
+ int dampflag = read_int(fp);
+
+ if (!flag) return;
+
+ for (i = 1; i <= data.ntypes; i++) {
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"lj/charmm/coul/charmm") == 0) ||
+ (strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) ||
+ (strcmp(style,"lj/charmm/coul/long") == 0)) {
+
+ if (strcmp(style,"lj/charmm/coul/charmm") == 0) {
+ double cut_lj_inner = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul_inner = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) {
+ double cut_lj_inner = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul_inner = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/charmm/coul/long") == 0) {
+ double cut_lj_inner = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ data.pair_charmm_epsilon = new double[data.ntypes+1];
+ data.pair_charmm_sigma = new double[data.ntypes+1];
+ data.pair_charmm_eps14 = new double[data.ntypes+1];
+ data.pair_charmm_sigma14 = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_charmm_epsilon[i] = read_double(fp);
+ data.pair_charmm_sigma[i] = read_double(fp);
+ data.pair_charmm_eps14[i] = read_double(fp);
+ data.pair_charmm_sigma14[i] = read_double(fp);
+ } else {
+ double charmm_epsilon = read_double(fp);
+ double charmm_sigma = read_double(fp);
+ double charmm_eps14 = read_double(fp);
+ double charmm_sigma14 = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"lj/class2") == 0) ||
+ (strcmp(style,"lj/class2/coul/cut") == 0) ||
+ (strcmp(style,"lj/class2/coul/long") == 0)) {
+
+ if (strcmp(style,"lj/class2") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/class2/coul/cut") == 0) {
+ m = 1;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/class2/coul/long") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ data.pair_class2_epsilon = new double[data.ntypes+1];
+ data.pair_class2_sigma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_class2_epsilon[i] = read_double(fp);
+ data.pair_class2_sigma[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ } else {
+ double class2_epsilon = read_double(fp);
+ double class2_sigma = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"lj/cut") == 0) ||
+ (strcmp(style,"lj96/cut") == 0) ||
+ (strcmp(style,"lj/cut/coul/cut") == 0) ||
+ (strcmp(style,"lj/cut/coul/debye") == 0) ||
+ (strcmp(style,"lj/cut/coul/long") == 0) ||
+ (strcmp(style,"lj/cut/coul/long/tip4p") == 0) ||
+ (strcmp(style,"lj/coul") == 0)) {
+
+ if (strcmp(style,"lj/cut") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj96/cut") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/cut/coul/cut") == 0) {
+ m = 1;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/cut/coul/debye") == 0) {
+ m = 1;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ double kappa = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/cut/coul/long") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/cut/coul/long/tip4p") == 0) {
+ m = 0;
+ int typeO = read_int(fp);
+ int typeH = read_int(fp);
+ int typeB = read_int(fp);
+ int typeA = read_int(fp);
+ double qdist = read_double(fp);
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/coul") == 0) {
+ m = 0;
+ double cut_lj_global = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ int ewald_order = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ data.pair_lj_epsilon = new double[data.ntypes+1];
+ data.pair_lj_sigma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_lj_epsilon[i] = read_double(fp);
+ data.pair_lj_sigma[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ } else {
+ double lj_epsilon = read_double(fp);
+ double lj_sigma = read_double(fp);
+ double cut_lj = read_double(fp);
+ if (m) double cut_coul = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"lj/expand") == 0) {
+
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_ljexpand_epsilon = new double[data.ntypes+1];
+ data.pair_ljexpand_sigma = new double[data.ntypes+1];
+ data.pair_ljexpand_shift = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_ljexpand_epsilon[i] = read_double(fp);
+ data.pair_ljexpand_sigma[i] = read_double(fp);
+ data.pair_ljexpand_shift[i] = read_double(fp);
+ double cut_lj = read_double(fp);
+ } else {
+ double ljexpand_epsilon = read_double(fp);
+ double ljexpand_sigma = read_double(fp);
+ double ljexpand_shift = read_double(fp);
+ double cut_lj = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"lj/gromacs") == 0) ||
+ (strcmp(style,"lj/gromacs/coul/gromacs") == 0)) {
+
+ if (strcmp(style,"lj/gromacs") == 0) {
+ m = 1;
+ double cut_inner_global = read_double(fp);
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ } else if (strcmp(style,"lj/gromacs/coul/gromacs") == 0) {
+ m = 0;
+ double cut_lj_inner_global = read_double(fp);
+ double cut_lj = read_double(fp);
+ double cut_coul_inner_global = read_double(fp);
+ double cut_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ }
+
+ if (!flag) return;
+
+ data.pair_ljgromacs_epsilon = new double[data.ntypes+1];
+ data.pair_ljgromacs_sigma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_ljgromacs_epsilon[i] = read_double(fp);
+ data.pair_ljgromacs_sigma[i] = read_double(fp);
+ if (m) {
+ double cut_inner = read_double(fp);
+ double cut = read_double(fp);
+ }
+ } else {
+ double ljgromacs_epsilon = read_double(fp);
+ double ljgromacs_sigma = read_double(fp);
+ if (m) {
+ double cut_inner = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+ }
+
+ } else if (strcmp(style,"lj/smooth") == 0) {
+
+ double cut_inner_global = read_double(fp);
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_ljsmooth_epsilon = new double[data.ntypes+1];
+ data.pair_ljsmooth_sigma = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_ljsmooth_epsilon[i] = read_double(fp);
+ data.pair_ljsmooth_sigma[i] = read_double(fp);
+ double cut_inner = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double ljsmooth_epsilon = read_double(fp);
+ double ljsmooth_sigma = read_double(fp);
+ double cut_inner = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"meam") == 0) {
+
+ } else if (strcmp(style,"morse") == 0) {
+
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_morse_d0 = new double[data.ntypes+1];
+ data.pair_morse_alpha = new double[data.ntypes+1];
+ data.pair_morse_r0 = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_morse_d0[i] = read_double(fp);
+ data.pair_morse_alpha[i] = read_double(fp);
+ data.pair_morse_r0[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double morse_d0 = read_double(fp);
+ double morse_alpha = read_double(fp);
+ double morse_r0 = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"reax") == 0) {
+ } else if (strcmp(style,"reax/c") == 0) {
+
+ } else if (strcmp(style,"soft") == 0) {
+
+ double cut_global = read_double(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_soft_A = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_soft_A[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double soft_A = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if (strcmp(style,"sw") == 0) {
+
+ } else if (strcmp(style,"table") == 0) {
+
+ int tabstyle = read_int(fp);
+ int n = read_int(fp);
+
+ } else if (strcmp(style,"tersoff") == 0) {
+ } else if (strcmp(style,"tersoff/zbl") == 0) {
+
+ } else if (strcmp(style,"yukawa") == 0) {
+
+ double kappa = read_double(fp);
+ double cut_global = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ data.pair_yukawa_A = new double[data.ntypes+1];
+
+ for (i = 1; i <= data.ntypes; i++)
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ if (i == j) {
+ data.pair_yukawa_A[i] = read_double(fp);
+ double cut = read_double(fp);
+ } else {
+ double yukawa_A = read_double(fp);
+ double cut = read_double(fp);
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"cg/cmm") == 0) ||
+ (strcmp(style,"cg/cmm/coul/cut") == 0) ||
+ (strcmp(style,"cg/cmm/coul/long") == 0)) {
+ m = 0;
+ data.cut_lj_global = read_double(fp);
+ data.cut_coul_global = read_double(fp);
+ data.kappa = read_double(fp);
+ data.offset_flag = read_int(fp);
+ data.mix_flag = read_int(fp);
+
+ if (!flag) return;
+
+ const int numtyp=data.ntypes+1;
+ data.pair_cg_cmm_type = new int*[numtyp];
+ data.pair_setflag = new int*[numtyp];
+ data.pair_cg_epsilon = new double*[numtyp];
+ data.pair_cg_sigma = new double*[numtyp];
+ data.pair_cut_lj = new double*[numtyp];
+ if ((strcmp(style,"cg/cmm/coul/cut") == 0) ||
+ (strcmp(style,"cg/cmm/coul/long") == 0)) {
+ data.pair_cut_coul = new double*[numtyp];
+ m=1;
+ } else {
+ data.pair_cut_coul = NULL;
+ m=0;
+ }
+
+ for (i = 1; i <= data.ntypes; i++) {
+ data.pair_cg_cmm_type[i] = new int[numtyp];
+ data.pair_setflag[i] = new int[numtyp];
+ data.pair_cg_epsilon[i] = new double[numtyp];
+ data.pair_cg_sigma[i] = new double[numtyp];
+ data.pair_cut_lj[i] = new double[numtyp];
+ if ((strcmp(style,"cg/cmm/coul/cut") == 0) ||
+ (strcmp(style,"cg/cmm/coul/long") == 0)) {
+ data.pair_cut_coul[i] = new double[numtyp];
+ }
+
+ for (j = i; j <= data.ntypes; j++) {
+ itmp = read_int(fp);
+ data.pair_setflag[i][j] = itmp;
+ if (i == j && itmp == 0) {
+ printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+ exit(1);
+ }
+ if (itmp) {
+ data.pair_cg_cmm_type[i][j] = read_int(fp);
+ data.pair_cg_epsilon[i][j] = read_double(fp);
+ data.pair_cg_sigma[i][j] = read_double(fp);
+ data.pair_cut_lj[i][j] = read_double(fp);
+ if (m) {
+ data.pair_cut_lj[i][j] = read_double(fp);
+ data.pair_cut_coul[i][j] = read_double(fp);
+ }
+ }
+ }
+ }
+
+ } else if ((strcmp(style,"hybrid") == 0) ||
+ (strcmp(style,"hybrid/overlay") == 0)) {
+
+ // for each substyle of hybrid,
+ // read its settings by calling pair() recursively with flag = 0
+ // so that coeff array allocation is skipped
+
+ int nstyles = read_int(fp);
+ for (int i = 0; i < nstyles; i++) {
+ char *substyle = read_char(fp);
+ pair(fp,data,substyle,0);
+ }
+
+ } else {
+ printf("ERROR: Unknown pair style %s\n",style);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// bond coeffs
+// one section for each bond style
+// ---------------------------------------------------------------------
+
+void bond(FILE *fp, Data &data)
+{
+ if (strcmp(data.bond_style,"none") == 0) {
+
+ } else if (strcmp(data.bond_style,"class2") == 0) {
+
+ data.bond_class2_r0 = new double[data.nbondtypes+1];
+ data.bond_class2_k2 = new double[data.nbondtypes+1];
+ data.bond_class2_k3 = new double[data.nbondtypes+1];
+ data.bond_class2_k4 = new double[data.nbondtypes+1];
+ fread(&data.bond_class2_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_class2_k2[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_class2_k3[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_class2_k4[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"fene") == 0) {
+
+ data.bond_fene_k = new double[data.nbondtypes+1];
+ data.bond_fene_r0 = new double[data.nbondtypes+1];
+ data.bond_fene_epsilon = new double[data.nbondtypes+1];
+ data.bond_fene_sigma = new double[data.nbondtypes+1];
+ fread(&data.bond_fene_k[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_fene_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_fene_epsilon[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_fene_sigma[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"fene/expand") == 0) {
+
+ data.bond_feneexpand_k = new double[data.nbondtypes+1];
+ data.bond_feneexpand_r0 = new double[data.nbondtypes+1];
+ data.bond_feneexpand_epsilon = new double[data.nbondtypes+1];
+ data.bond_feneexpand_sigma = new double[data.nbondtypes+1];
+ data.bond_feneexpand_shift = new double[data.nbondtypes+1];
+ fread(&data.bond_feneexpand_k[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_feneexpand_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_feneexpand_epsilon[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_feneexpand_sigma[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_feneexpand_shift[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"harmonic") == 0) {
+
+ data.bond_harmonic_k = new double[data.nbondtypes+1];
+ data.bond_harmonic_r0 = new double[data.nbondtypes+1];
+ fread(&data.bond_harmonic_k[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_harmonic_r0[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"harmonicshift") == 0) {
+
+ data.bond_harmonicshift_umin = new double[data.nbondtypes+1];
+ data.bond_harmonicshift_r0 = new double[data.nbondtypes+1];
+ data.bond_harmonicshift_rc = new double[data.nbondtypes+1];
+ fread(&data.bond_harmonicshift_umin[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_harmonicshift_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_harmonicshift_rc[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"harmonicshiftcut") == 0) {
+
+ data.bond_harmonicshiftcut_umin = new double[data.nbondtypes+1];
+ data.bond_harmonicshiftcut_r0 = new double[data.nbondtypes+1];
+ data.bond_harmonicshiftcut_rc = new double[data.nbondtypes+1];
+ fread(&data.bond_harmonicshiftcut_umin[1],sizeof(double),
+ data.nbondtypes,fp);
+ fread(&data.bond_harmonicshiftcut_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_harmonicshiftcut_rc[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"morse") == 0) {
+
+ data.bond_morse_d0 = new double[data.nbondtypes+1];
+ data.bond_morse_alpha = new double[data.nbondtypes+1];
+ data.bond_morse_r0 = new double[data.nbondtypes+1];
+ fread(&data.bond_morse_d0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_morse_alpha[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_morse_r0[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"nonlinear") == 0) {
+
+ data.bond_nonlinear_epsilon = new double[data.nbondtypes+1];
+ data.bond_nonlinear_r0 = new double[data.nbondtypes+1];
+ data.bond_nonlinear_lamda = new double[data.nbondtypes+1];
+ fread(&data.bond_nonlinear_epsilon[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_nonlinear_r0[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_nonlinear_lamda[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"quartic") == 0) {
+
+ data.bond_quartic_k = new double[data.nbondtypes+1];
+ data.bond_quartic_b1 = new double[data.nbondtypes+1];
+ data.bond_quartic_b2 = new double[data.nbondtypes+1];
+ data.bond_quartic_rc = new double[data.nbondtypes+1];
+ data.bond_quartic_u0 = new double[data.nbondtypes+1];
+ fread(&data.bond_quartic_k[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_quartic_b1[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_quartic_b2[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_quartic_rc[1],sizeof(double),data.nbondtypes,fp);
+ fread(&data.bond_quartic_u0[1],sizeof(double),data.nbondtypes,fp);
+
+ } else if (strcmp(data.bond_style,"table") == 0) {
+
+ int tabstyle = read_int(fp);
+ int n = read_int(fp);
+
+ } else if (strcmp(data.bond_style,"hybrid") == 0) {
+
+ int nstyles = read_int(fp);
+ for (int i = 0; i < nstyles; i++)
+ char *substyle = read_char(fp);
+
+ } else {
+ printf("ERROR: Unknown bond style %s\n",data.bond_style);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// angle coeffs
+// one section for each angle style
+// ---------------------------------------------------------------------
+
+void angle(FILE *fp, Data &data)
+{
+ if (strcmp(data.angle_style,"none") == 0) {
+
+ } else if (strcmp(data.angle_style,"cg/cmm") == 0) {
+
+ data.angle_harmonic_k = new double[data.nangletypes+1];
+ data.angle_harmonic_theta0 = new double[data.nangletypes+1];
+ data.angle_cg_cmm_epsilon = new double[data.nangletypes+1];
+ data.angle_cg_cmm_sigma = new double[data.nangletypes+1];
+ double *angle_cg_cmm_rcut = new double[data.nangletypes+1];
+ data.angle_cg_cmm_type = new int[data.nangletypes+1];
+
+ fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_cg_cmm_epsilon[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_cg_cmm_sigma[1],sizeof(double),data.nangletypes,fp);
+ fread(angle_cg_cmm_rcut,sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_cg_cmm_type[1],sizeof(int),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"charmm") == 0) {
+
+ data.angle_charmm_k = new double[data.nangletypes+1];
+ data.angle_charmm_theta0 = new double[data.nangletypes+1];
+ data.angle_charmm_k_ub = new double[data.nangletypes+1];
+ data.angle_charmm_r_ub = new double[data.nangletypes+1];
+ fread(&data.angle_charmm_k[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_charmm_theta0[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_charmm_k_ub[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_charmm_r_ub[1],sizeof(double),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"class2") == 0) {
+
+ data.angle_class2_theta0 = new double[data.nangletypes+1];
+ data.angle_class2_k2 = new double[data.nangletypes+1];
+ data.angle_class2_k3 = new double[data.nangletypes+1];
+ data.angle_class2_k4 = new double[data.nangletypes+1];
+
+ data.angle_class2_bb_k = new double[data.nangletypes+1];
+ data.angle_class2_bb_r1 = new double[data.nangletypes+1];
+ data.angle_class2_bb_r2 = new double[data.nangletypes+1];
+
+ data.angle_class2_ba_k1 = new double[data.nangletypes+1];
+ data.angle_class2_ba_k2 = new double[data.nangletypes+1];
+ data.angle_class2_ba_r1 = new double[data.nangletypes+1];
+ data.angle_class2_ba_r2 = new double[data.nangletypes+1];
+
+ fread(&data.angle_class2_theta0[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_k2[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_k3[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_k4[1],sizeof(double),data.nangletypes,fp);
+
+ fread(&data.angle_class2_bb_k[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_bb_r1[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_bb_r2[1],sizeof(double),data.nangletypes,fp);
+
+ fread(&data.angle_class2_ba_k1[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_ba_k2[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_ba_r1[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_class2_ba_r2[1],sizeof(double),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"cosineshift") == 0) {
+
+ data.angle_cosineshift_umin = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshift_umin[1],sizeof(double),data.nangletypes,fp);
+ data.angle_cosineshift_cost = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshift_cost[1],sizeof(double),data.nangletypes,fp);
+ data.angle_cosineshift_sint = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshift_sint[1],sizeof(double),data.nangletypes,fp);
+ data.angle_cosineshift_theta0 = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshift_theta0[1],sizeof(double),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"cosineshiftexp") == 0) {
+
+ data.angle_cosineshiftexp_umin = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshiftexp_umin[1],sizeof(double),
+ data.nangletypes,fp);
+ data.angle_cosineshiftexp_a = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshiftexp_a[1],sizeof(double),data.nangletypes,fp);
+ data.angle_cosineshiftexp_cost = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshiftexp_cost[1],sizeof(double),
+ data.nangletypes,fp);
+ data.angle_cosineshiftexp_sint = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshiftexp_sint[1],sizeof(double),
+ data.nangletypes,fp);
+ data.angle_cosineshiftexp_theta0 = new double[data.nangletypes+1];
+ fread(&data.angle_cosineshiftexp_theta0[1],sizeof(double),
+ data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"cosine") == 0) {
+
+ data.angle_cosine_k = new double[data.nangletypes+1];
+ fread(&data.angle_cosine_k[1],sizeof(double),data.nangletypes,fp);
+
+ } else if ((strcmp(data.angle_style,"cosine/squared") == 0) ||
+ (strcmp(data.angle_style,"cosine/delta") == 0)) {
+
+ data.angle_cosine_squared_k = new double[data.nangletypes+1];
+ data.angle_cosine_squared_theta0 = new double[data.nangletypes+1];
+ fread(&data.angle_cosine_squared_k[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_cosine_squared_theta0[1],
+ sizeof(double),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"harmonic") == 0) {
+
+ data.angle_harmonic_k = new double[data.nangletypes+1];
+ data.angle_harmonic_theta0 = new double[data.nangletypes+1];
+ fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp);
+ fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp);
+
+ } else if (strcmp(data.angle_style,"table") == 0) {
+
+ int tabstyle = read_int(fp);
+ int n = read_int(fp);
+
+ } else if (strcmp(data.angle_style,"hybrid") == 0) {
+
+ int nstyles = read_int(fp);
+ for (int i = 0; i < nstyles; i++)
+ char *substyle = read_char(fp);
+
+ } else {
+ printf("ERROR: Unknown angle style %s\n",data.angle_style);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// dihedral coeffs
+// one section for each dihedral style
+// ---------------------------------------------------------------------
+
+void dihedral(FILE *fp, Data &data)
+{
+ if (strcmp(data.dihedral_style,"none") == 0) {
+
+ } else if (strcmp(data.dihedral_style,"charmm") == 0) {
+
+ data.dihedral_charmm_k = new double[data.ndihedraltypes+1];
+ data.dihedral_charmm_multiplicity = new int[data.ndihedraltypes+1];
+ data.dihedral_charmm_sign = new int[data.ndihedraltypes+1];
+ data.dihedral_charmm_weight = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_charmm_k[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_charmm_multiplicity[1],sizeof(int),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_charmm_sign[1],sizeof(int),data.ndihedraltypes,fp);
+ fread(&data.dihedral_charmm_weight[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"class2") == 0) {
+
+ data.dihedral_class2_k1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_k2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_k3 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_phi1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_phi2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_phi3 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_mbt_f1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_mbt_f2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_mbt_f3 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_mbt_r0 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_ebt_f1_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_f2_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_f3_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_r0_1 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_ebt_f1_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_f2_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_f3_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_ebt_r0_2 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_at_f1_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_f2_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_f3_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_theta0_1 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_at_f1_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_f2_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_f3_2 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_at_theta0_2 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_aat_k = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_aat_theta0_1 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_aat_theta0_2 = new double[data.ndihedraltypes+1];
+
+ data.dihedral_class2_bb13_k = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_bb13_r10 = new double[data.ndihedraltypes+1];
+ data.dihedral_class2_bb13_r30 = new double[data.ndihedraltypes+1];
+
+ fread(&data.dihedral_class2_k1[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_k2[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_k3[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_phi1[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_phi2[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_phi3[1],sizeof(double),data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_mbt_f1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_mbt_f2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_mbt_f3[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_mbt_r0[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_ebt_f1_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_f2_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_f3_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_r0_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_ebt_f1_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_f2_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_f3_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_ebt_r0_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_at_f1_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_f2_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_f3_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_theta0_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_at_f1_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_f2_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_f3_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_at_theta0_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_aat_k[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_aat_theta0_1[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_aat_theta0_2[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ fread(&data.dihedral_class2_bb13_k[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_bb13_r10[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_class2_bb13_r30[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"harmonic") == 0) {
+
+ data.dihedral_harmonic_k = new double[data.ndihedraltypes+1];
+ data.dihedral_harmonic_multiplicity = new int[data.ndihedraltypes+1];
+ data.dihedral_harmonic_sign = new int[data.ndihedraltypes+1];
+ fread(&data.dihedral_harmonic_k[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_harmonic_multiplicity[1],sizeof(int),
+ data.ndihedraltypes,fp);
+ fread(&data.dihedral_harmonic_sign[1],sizeof(int),data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"helix") == 0) {
+
+ data.dihedral_helix_aphi = new double[data.ndihedraltypes+1];
+ data.dihedral_helix_bphi = new double[data.ndihedraltypes+1];
+ data.dihedral_helix_cphi = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_helix_aphi[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_helix_bphi[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_helix_cphi[1],sizeof(double),data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"multi/harmonic") == 0) {
+
+ data.dihedral_multi_a1 = new double[data.ndihedraltypes+1];
+ data.dihedral_multi_a2 = new double[data.ndihedraltypes+1];
+ data.dihedral_multi_a3 = new double[data.ndihedraltypes+1];
+ data.dihedral_multi_a4 = new double[data.ndihedraltypes+1];
+ data.dihedral_multi_a5 = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_multi_a1[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_multi_a2[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_multi_a3[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_multi_a4[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_multi_a5[1],sizeof(double),data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"opls") == 0) {
+
+ data.dihedral_opls_k1 = new double[data.ndihedraltypes+1];
+ data.dihedral_opls_k2 = new double[data.ndihedraltypes+1];
+ data.dihedral_opls_k3 = new double[data.ndihedraltypes+1];
+ data.dihedral_opls_k4 = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_opls_k1[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_opls_k2[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_opls_k3[1],sizeof(double),data.ndihedraltypes,fp);
+ fread(&data.dihedral_opls_k4[1],sizeof(double),data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"cosineshiftexp") == 0) {
+
+ data.dihedral_cosineshiftexp_umin = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_cosineshiftexp_umin[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ data.dihedral_cosineshiftexp_a = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_cosineshiftexp_a[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ data.dihedral_cosineshiftexp_cost = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_cosineshiftexp_cost[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ data.dihedral_cosineshiftexp_sint = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_cosineshiftexp_sint[1],sizeof(double),
+ data.ndihedraltypes,fp);
+ data.dihedral_cosineshiftexp_theta = new double[data.ndihedraltypes+1];
+ fread(&data.dihedral_cosineshiftexp_theta[1],sizeof(double),
+ data.ndihedraltypes,fp);
+
+ } else if (strcmp(data.dihedral_style,"table") == 0) {
+
+ int tabstyle = read_int(fp);
+ int n = read_int(fp);
+
+ } else if (strcmp(data.dihedral_style,"hybrid") == 0) {
+
+ int nstyles = read_int(fp);
+ for (int i = 0; i < nstyles; i++)
+ char *substyle = read_char(fp);
+
+ } else {
+ printf("ERROR: Unknown dihedral style %s\n",data.dihedral_style);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// improper coeffs
+// one section for each improper style
+// ---------------------------------------------------------------------
+
+void improper(FILE *fp, Data &data)
+{
+ if (strcmp(data.improper_style,"none") == 0) {
+
+ } else if (strcmp(data.improper_style,"class2") == 0) {
+
+ data.improper_class2_k0 = new double[data.nimpropertypes+1];
+ data.improper_class2_chi0 = new double[data.nimpropertypes+1];
+
+ data.improper_class2_aa_k1 = new double[data.nimpropertypes+1];
+ data.improper_class2_aa_k2 = new double[data.nimpropertypes+1];
+ data.improper_class2_aa_k3 = new double[data.nimpropertypes+1];
+ data.improper_class2_aa_theta0_1 = new double[data.nimpropertypes+1];
+ data.improper_class2_aa_theta0_2 = new double[data.nimpropertypes+1];
+ data.improper_class2_aa_theta0_3 = new double[data.nimpropertypes+1];
+
+ fread(&data.improper_class2_k0[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_chi0[1],sizeof(double),
+ data.nimpropertypes,fp);
+
+ fread(&data.improper_class2_aa_k1[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_aa_k2[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_aa_k3[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_aa_theta0_1[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_aa_theta0_2[1],sizeof(double),
+ data.nimpropertypes,fp);
+ fread(&data.improper_class2_aa_theta0_3[1],sizeof(double),
+ data.nimpropertypes,fp);
+
+ } else if (strcmp(data.improper_style,"cvff") == 0) {
+
+ data.improper_cvff_k = new double[data.nimpropertypes+1];
+ data.improper_cvff_sign = new int[data.nimpropertypes+1];
+ data.improper_cvff_multiplicity = new int[data.nimpropertypes+1];
+ fread(&data.improper_cvff_k[1],sizeof(double),data.nimpropertypes,fp);
+ fread(&data.improper_cvff_sign[1],sizeof(int),data.nimpropertypes,fp);
+ fread(&data.improper_cvff_multiplicity[1],sizeof(int),
+ data.nimpropertypes,fp);
+
+ } else if (strcmp(data.improper_style,"harmonic") == 0) {
+
+ data.improper_harmonic_k = new double[data.nimpropertypes+1];
+ data.improper_harmonic_chi = new double[data.nimpropertypes+1];
+ fread(&data.improper_harmonic_k[1],sizeof(double),data.nimpropertypes,fp);
+ fread(&data.improper_harmonic_chi[1],sizeof(double),
+ data.nimpropertypes,fp);
+
+ } else if (strcmp(data.improper_style,"hybrid") == 0) {
+
+ int nstyles = read_int(fp);
+ for (int i = 0; i < nstyles; i++)
+ char *substyle = read_char(fp);
+
+ } else {
+ printf("ERROR: Unknown improper style %s\n",data.improper_style);
+ exit(1);
+ }
+}
+
+// ---------------------------------------------------------------------
+// initialize Data
+// ---------------------------------------------------------------------
+
+Data::Data()
+{
+ nellipsoids = 0;
+}
+
+// ---------------------------------------------------------------------
+// print out stats on problem
+// ---------------------------------------------------------------------
+
+void Data::stats()
+{
+ char fstr[64];
+
+ printf(" Restart file version = %s\n",version);
+ printf(" Ntimestep = " BIGINT_FORMAT "\n",ntimestep);
+ printf(" Nprocs = %d\n",nprocs);
+ printf(" Natoms = " BIGINT_FORMAT "\n",natoms);
+
+ if (nellipsoids) printf(" Nellipsoids = " BIGINT_FORMAT "\n",nellipsoids);
+
+ printf(" Nbonds = " BIGINT_FORMAT "\n",nbonds);
+ printf(" Nangles = " BIGINT_FORMAT "\n",nangles);
+ printf(" Ndihedrals = " BIGINT_FORMAT "\n",ndihedrals);
+ printf(" Nimpropers = " BIGINT_FORMAT "\n",nimpropers);
+
+ printf(" Unit style = %s\n",unit_style);
+ printf(" Atom style = %s\n",atom_style);
+ printf(" Pair style = %s\n",pair_style);
+ printf(" Bond style = %s\n",bond_style);
+ printf(" Angle style = %s\n",angle_style);
+ printf(" Dihedral style = %s\n",dihedral_style);
+ printf(" Improper style = %s\n",improper_style);
+
+ printf(" Xlo xhi = %g %g\n",xlo,xhi);
+ printf(" Ylo yhi = %g %g\n",ylo,yhi);
+ printf(" Zlo zhi = %g %g\n",zlo,zhi);
+ if (triclinic) printf(" Xy xz yz = %g %g %g\n",xy,xz,yz);
+ printf(" Periodicity = %d %d %d\n",xperiodic,yperiodic,zperiodic);
+ printf(" Boundary = %d %d, %d %d, %d %d\n",boundary[0][0],boundary[0][1],
+ boundary[1][0],boundary[1][1],boundary[2][0],boundary[2][1]);
+}
+
+// ---------------------------------------------------------------------
+// write the data file and input file
+// ---------------------------------------------------------------------
+
+void Data::write(FILE *fp, FILE *fp2)
+{
+ fprintf(fp,"LAMMPS data file from restart file: timestep = "
+ BIGINT_FORMAT ", procs = %d\n\n",ntimestep,nprocs);
+
+ fprintf(fp,BIGINT_FORMAT " atoms\n",natoms);
+ if (nellipsoids) fprintf(fp,BIGINT_FORMAT " ellipsoids\n",nellipsoids);
+ if (nbonds) fprintf(fp,BIGINT_FORMAT " bonds\n",nbonds);
+ if (nangles) fprintf(fp,BIGINT_FORMAT " angles\n",nangles);
+ if (ndihedrals) fprintf(fp,BIGINT_FORMAT " dihedrals\n",ndihedrals);
+ if (nimpropers) fprintf(fp,BIGINT_FORMAT " impropers\n",nimpropers);
+
+ fprintf(fp,"\n");
+
+ fprintf(fp,"%d atom types\n",ntypes);
+ if (nbondtypes) fprintf(fp,"%d bond types\n",nbondtypes);
+ if (nangletypes) fprintf(fp,"%d angle types\n",nangletypes);
+ if (ndihedraltypes) fprintf(fp,"%d dihedral types\n",ndihedraltypes);
+ if (nimpropertypes) fprintf(fp,"%d improper types\n",nimpropertypes);
+
+ fprintf(fp,"\n");
+
+ fprintf(fp,"%-1.16e %-1.16e xlo xhi\n",xlo,xhi);
+ fprintf(fp,"%-1.16e %-1.16e ylo yhi\n",ylo,yhi);
+ fprintf(fp,"%-1.16e %-1.16e zlo zhi\n",zlo,zhi);
+ if (triclinic) fprintf(fp,"%-1.16e %-1.16e %-1.16e xy xz yz\n",xy,xz,yz);
+
+ // write ff styles to input file
+
+ if (fp2) {
+ fprintf(fp2,"# LAMMPS input file from restart file: "
+ "timestep = " BIGINT_FORMAT ", procs = %d\n\n",ntimestep,nprocs);
+
+ if (pair_style) fprintf(fp2,"pair_style %s\n",pair_style);
+ if (bond_style) fprintf(fp2,"bond_style %s\n",bond_style);
+ if (angle_style) fprintf(fp2,"angle_style %s\n",angle_style);
+ if (dihedral_style) fprintf(fp2,"dihedral_style %s\n",dihedral_style);
+ if (improper_style) fprintf(fp2,"improper_style %s\n",improper_style);
+ fprintf(fp2,"special_bonds %g %g %g %g %g %g\n",
+ special_lj[1],special_lj[2],special_lj[3],
+ special_lj[1],special_coul[2],special_coul[3]);
+ fprintf(fp2,"\n");
+ }
+
+ // mass to either data file or input file
+
+ if (mass) {
+ if (fp2) {
+ fprintf(fp2,"\n");
+ for (int i = 1; i <= ntypes; i++) fprintf(fp2,"mass %d %g\n",i,mass[i]);
+ fprintf(fp2,"\n");
+ } else {
+ fprintf(fp,"\nMasses\n\n");
+ for (int i = 1; i <= ntypes; i++) fprintf(fp,"%d %g\n",i,mass[i]);
+ }
+ }
+
+ // pair coeffs to data file
+
+ if (pair_style && fp2 == NULL) {
+ if ((strcmp(pair_style,"none") != 0) &&
+ (strcmp(pair_style,"adp") != 0) &&
+ (strcmp(pair_style,"airebo") != 0) &&
+ (strcmp(pair_style,"coul/cut") != 0) &&
+ (strcmp(pair_style,"coul/debye") != 0) &&
+ (strcmp(pair_style,"coul/long") != 0) &&
+ (strcmp(pair_style,"eam") != 0) &&
+ (strcmp(pair_style,"eam/alloy") != 0) &&
+ (strcmp(pair_style,"eam/fs") != 0) &&
+ (strcmp(pair_style,"eim") != 0) &&
+ (strcmp(pair_style,"eff/cut") != 0) &&
+ (strcmp(pair_style,"gran/history") != 0) &&
+ (strcmp(pair_style,"gran/no_history") != 0) &&
+ (strcmp(pair_style,"gran/hertzian") != 0) &&
+ (strcmp(pair_style,"meam") != 0) &&
+ (strcmp(pair_style,"reax") != 0) &&
+ (strcmp(pair_style,"reax/c") != 0) &&
+ (strcmp(pair_style,"sw") != 0) &&
+ (strcmp(pair_style,"table") != 0) &&
+ (strcmp(pair_style,"tersoff") != 0) &&
+ (strcmp(pair_style,"tersoff/zbl") != 0) &&
+ (strcmp(pair_style,"hybrid") != 0) &&
+ (strcmp(pair_style,"hybrid/overlay") != 0))
+ fprintf(fp,"\nPair Coeffs\n\n");
+
+ if (strcmp(pair_style,"born/coul/long") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g\n",i,
+ pair_born_A[i],pair_born_rho[i],pair_born_sigma[i],
+ pair_born_C[i],pair_born_D[i]);
+
+ } else if ((strcmp(pair_style,"buck") == 0) ||
+ (strcmp(pair_style,"buck/coul/cut") == 0) ||
+ (strcmp(pair_style,"buck/coul/long") == 0) ||
+ (strcmp(pair_style,"buck/long") == 0)) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ pair_buck_A[i],pair_buck_rho[i],pair_buck_C[i]);
+
+ } else if (strcmp(pair_style,"colloid") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ pair_colloid_A12[i],pair_colloid_sigma[i],
+ pair_colloid_d2[i],pair_colloid_d2[i]);
+
+ } else if (strcmp(pair_style,"dipole/cut") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_dipole_epsilon[i],pair_dipole_sigma[i]);
+
+ } else if (strcmp(pair_style,"dpd") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_dpd_a0[i],pair_dpd_gamma[i]);
+
+ } else if (strcmp(pair_style,"dpd/tstat") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g\n",i,
+ pair_dpd_gamma[i]);
+
+ } else if (strcmp(pair_style,"gayberne") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+ pair_gb_epsilon[i],pair_gb_sigma[i],
+ pair_gb_epsa[i],pair_gb_epsb[i],pair_gb_epsc[i],
+ pair_gb_epsa[i],pair_gb_epsb[i],pair_gb_epsc[i]);
+
+ } else if ((strcmp(pair_style,"lj/charmm/coul/charmm") == 0) ||
+ (strcmp(pair_style,"lj/charmm/coul/charmm/implicit") == 0) ||
+ (strcmp(pair_style,"lj/charmm/coul/long") == 0)) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ pair_charmm_epsilon[i],pair_charmm_sigma[i],
+ pair_charmm_eps14[i],pair_charmm_sigma14[i]);
+
+ } else if ((strcmp(pair_style,"lj/class2") == 0) ||
+ (strcmp(pair_style,"lj/class2/coul/cut") == 0) ||
+ (strcmp(pair_style,"lj/class2/coul/long") == 0)) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_class2_epsilon[i],pair_class2_sigma[i]);
+
+ } else if ((strcmp(pair_style,"lj/cut") == 0) ||
+ (strcmp(pair_style,"lj96/cut") == 0) ||
+ (strcmp(pair_style,"lj/cut/coul/cut") == 0) ||
+ (strcmp(pair_style,"lj/cut/coul/debye") == 0) ||
+ (strcmp(pair_style,"lj/cut/coul/long") == 0) ||
+ (strcmp(pair_style,"lj/cut/coul/long/tip4p") == 0) ||
+ (strcmp(pair_style,"lj/coul") == 0)) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_lj_epsilon[i],pair_lj_sigma[i]);
+
+ } else if (strcmp(pair_style,"lj/expand") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ pair_ljexpand_epsilon[i],pair_ljexpand_sigma[i],
+ pair_ljexpand_shift[i]);
+
+ } else if ((strcmp(pair_style,"lj/gromacs") == 0) ||
+ (strcmp(pair_style,"lj/gromacs/coul/gromacs") == 0)) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_ljgromacs_epsilon[i],pair_ljgromacs_sigma[i]);
+
+ } else if (strcmp(pair_style,"lj/smooth") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ pair_ljsmooth_epsilon[i],pair_ljsmooth_sigma[i]);
+
+ } else if (strcmp(pair_style,"morse") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ pair_morse_d0[i],pair_morse_alpha[i],pair_morse_r0[i]);
+
+ } else if (strcmp(pair_style,"soft") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g\n",i,
+ pair_soft_A[i]);
+
+ } else if (strcmp(pair_style,"yukawa") == 0) {
+ for (int i = 1; i <= ntypes; i++)
+ fprintf(fp,"%d %g\n",i,
+ pair_yukawa_A[i]);
+
+ } else if ((strcmp(pair_style,"cg/cmm") == 0) ||
+ (strcmp(pair_style,"cg/cmm/coul/cut") == 0) ||
+ (strcmp(pair_style,"cg/cmm/coul/long") == 0)) {
+ printf("ERROR: Cannot write pair_style %s to data file\n",
+ pair_style);
+ exit(1);
+ }
+ }
+
+ // pair coeffs to input file
+ // only supported styles = cg/cmm
+
+ if (pair_style && fp2) {
+ if ((strcmp(pair_style,"cg/cmm") == 0) ||
+ (strcmp(pair_style,"cg/cmm/coul/cut") == 0) ||
+ (strcmp(pair_style,"cg/cmm/coul/long") == 0)) {
+ for (int i = 1; i <= ntypes; i++) {
+ for (int j = i; j <= ntypes; j++) {
+ fprintf(fp2,"pair_coeff %d %d %s %g %g\n",i,j,
+ cg_type_list[pair_cg_cmm_type[i][j]],
+ pair_cg_epsilon[i][j],pair_cg_sigma[i][j]);
+ }
+ }
+
+ } else {
+ printf("ERROR: Cannot write pair_style %s to input file\n",
+ pair_style);
+ exit(1);
+ }
+ }
+
+ // bond coeffs to data file
+
+ if (bond_style && fp2 == NULL) {
+ if ((strcmp(bond_style,"none") != 0) &&
+ (strcmp(bond_style,"table") != 0) &&
+ (strcmp(bond_style,"hybrid") != 0))
+ fprintf(fp,"\nBond Coeffs\n\n");
+
+ if (strcmp(bond_style,"class2") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ bond_class2_r0[i],bond_class2_k2[i],
+ bond_class2_k3[i],bond_class2_k4[i]);
+
+ } else if (strcmp(bond_style,"fene") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ bond_fene_k[i],bond_fene_r0[i],
+ bond_fene_epsilon[i],bond_fene_sigma[i]);
+
+ } else if (strcmp(bond_style,"fene/expand") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g\n",i,
+ bond_feneexpand_k[i],bond_feneexpand_r0[i],
+ bond_feneexpand_epsilon[i],bond_feneexpand_sigma[i],
+ bond_feneexpand_shift[i]);
+
+ } else if (strcmp(bond_style,"harmonic") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ bond_harmonic_k[i],bond_harmonic_r0[i]);
+
+ } else if (strcmp(bond_style,"harmonicshift") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ bond_harmonicshift_umin[i],bond_harmonicshift_r0[i],
+ bond_harmonicshift_rc[i]);
+
+ } else if (strcmp(bond_style,"harmonicshiftcut") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ bond_harmonicshiftcut_umin[i],bond_harmonicshiftcut_r0[i],
+ bond_harmonicshiftcut_rc[i]);
+
+ } else if (strcmp(bond_style,"morse") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ bond_morse_d0[i],bond_morse_alpha[i],bond_morse_r0[i]);
+
+ } else if (strcmp(bond_style,"nonlinear") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ bond_nonlinear_epsilon[i],bond_nonlinear_r0[i],
+ bond_nonlinear_lamda[i]);
+
+ } else if (strcmp(bond_style,"quartic") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g\n",i,
+ bond_quartic_k[i],bond_quartic_b1[i],bond_quartic_b2[i],
+ bond_quartic_rc[i],bond_quartic_u0[i]);
+ }
+ }
+
+ // bond coeffs to input file
+ // only supported styles = harmonic
+
+ if (bond_style && fp2) {
+ if (strcmp(bond_style,"harmonic") == 0) {
+ for (int i = 1; i <= nbondtypes; i++)
+ fprintf(fp2,"bond_coeff %d %g %g\n",i,
+ bond_harmonic_k[i],bond_harmonic_r0[i]);
+
+ } else {
+ printf("ERROR: Cannot write bond_style %s to input file\n",
+ bond_style);
+ exit(1);
+ }
+ }
+
+ // angle coeffs to data file
+
+ if (angle_style && fp2 == NULL) {
+ if ((strcmp(angle_style,"none") != 0) &&
+ (strcmp(angle_style,"table") != 0) &&
+ (strcmp(angle_style,"hybrid") != 0))
+ fprintf(fp,"\nAngle Coeffs\n\n");
+
+ if (strcmp(angle_style,"charmm") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ angle_charmm_k[i],angle_charmm_theta0[i]/PI*180.0,
+ angle_charmm_k_ub[i],angle_charmm_r_ub[i]);
+
+ } else if (strcmp(angle_style,"class2") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ angle_class2_theta0[i]/PI*180.0,angle_class2_k2[i],
+ angle_class2_k3[i],angle_class2_k4[i]);
+
+ fprintf(fp,"\nBondBond Coeffs\n\n");
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ angle_class2_bb_k[i],
+ angle_class2_bb_r1[i],angle_class2_bb_r2[i]);
+
+ fprintf(fp,"\nBondAngle Coeffs\n\n");
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ angle_class2_ba_k1[i],angle_class2_ba_k2[i],
+ angle_class2_ba_r1[i],angle_class2_ba_r2[i]);
+
+ } else if (strcmp(angle_style,"cosine") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g\n",i,angle_cosine_k[i]);
+
+ } else if (strcmp(angle_style,"cosineshift") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g\n",i,angle_cosineshift_umin[i],
+ angle_cosineshift_theta0[i]/PI*180.0);
+
+ } else if (strcmp(angle_style,"cosineshiftexp") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,angle_cosineshiftexp_umin[i],
+ angle_cosineshiftexp_theta0[i]/PI*180.0,
+ angle_cosineshiftexp_a[i]);
+
+ } else if ((strcmp(angle_style,"cosine/squared") == 0) ||
+ (strcmp(angle_style,"cosine/delta") == 0)) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ angle_cosine_squared_k[i],
+ angle_cosine_squared_theta0[i]/PI*180.0);
+
+ } else if (strcmp(angle_style,"harmonic") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0);
+
+ } else if (strcmp(angle_style,"cg/cmm") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp,"%d %g %g %s %g %g\n",i,
+ angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0,
+ cg_type_list[angle_cg_cmm_type[i]],angle_cg_cmm_epsilon[i],
+ angle_cg_cmm_sigma[i]);
+ }
+ }
+
+ // angle coeffs to input file
+ // only supported styles = cosine/squared, harmonic, cg/cmm
+
+ if (angle_style && fp2) {
+ if ((strcmp(angle_style,"cosine/squared") == 0) ||
+ (strcmp(angle_style,"cosine/delta") == 0)) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp2,"angle_coeffs %d %g %g\n",i,
+ angle_cosine_squared_k[i],
+ angle_cosine_squared_theta0[i]/PI*180.0);
+
+ } else if (strcmp(angle_style,"harmonic") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp2,"angle_coeffs %d %g %g\n",i,
+ angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0);
+
+ } else if (strcmp(angle_style,"cg/cmm") == 0) {
+ for (int i = 1; i <= nangletypes; i++)
+ fprintf(fp2,"angle_coeffs %d %g %g %s %g %g\n",i,
+ angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0,
+ cg_type_list[angle_cg_cmm_type[i]],angle_cg_cmm_epsilon[i],
+ angle_cg_cmm_sigma[i]);
+
+ } else {
+ printf("ERROR: Cannot write angle_style %s to input file\n",
+ angle_style);
+ exit(1);
+ }
+ }
+
+ if (dihedral_style) {
+ if ((strcmp(dihedral_style,"none") != 0) &&
+ (strcmp(dihedral_style,"table") != 0) &&
+ (strcmp(dihedral_style,"hybrid") != 0))
+ fprintf(fp,"\nDihedral Coeffs\n\n");
+
+ if (strcmp(dihedral_style,"charmm") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %d %d %g\n",i,
+ dihedral_charmm_k[i],dihedral_charmm_multiplicity[i],
+ dihedral_charmm_sign[i],dihedral_charmm_weight[i]);
+
+ } else if (strcmp(dihedral_style,"class2") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g\n",i,
+ dihedral_class2_k1[i],
+ dihedral_class2_phi1[i]/PI*180.0,
+ dihedral_class2_k2[i],
+ dihedral_class2_phi2[i]/PI*180.0,
+ dihedral_class2_k3[i],
+ dihedral_class2_phi3[i]/PI*180.0);
+
+ fprintf(fp,"\nMiddleBondTorsion Coeffs\n\n");
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ dihedral_class2_mbt_f1[i],dihedral_class2_mbt_f2[i],
+ dihedral_class2_mbt_f3[i],dihedral_class2_mbt_r0[i]);
+
+ fprintf(fp,"\nEndBondTorsion Coeffs\n\n");
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+ dihedral_class2_ebt_f1_1[i],dihedral_class2_ebt_f2_1[i],
+ dihedral_class2_ebt_f3_1[i],
+ dihedral_class2_ebt_f1_2[i],dihedral_class2_ebt_f2_2[i],
+ dihedral_class2_ebt_f3_2[i],
+ dihedral_class2_ebt_r0_1[i],
+ dihedral_class2_ebt_r0_2[i]);
+
+ fprintf(fp,"\nAngleTorsion Coeffs\n\n");
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+ dihedral_class2_at_f1_1[i],dihedral_class2_at_f2_1[i],
+ dihedral_class2_at_f3_1[i],
+ dihedral_class2_at_f1_2[i],dihedral_class2_at_f2_2[i],
+ dihedral_class2_at_f3_2[i],
+ dihedral_class2_at_theta0_1[i]/PI*180.0,
+ dihedral_class2_at_theta0_2[i]/PI*180.0);
+
+ fprintf(fp,"\nAngleAngleTorsion Coeffs\n\n");
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ dihedral_class2_aat_k[i],
+ dihedral_class2_aat_theta0_1[i]/PI*180.0,
+ dihedral_class2_aat_theta0_2[i]/PI*180.0);
+
+ fprintf(fp,"\nBondBond13 Coeffs\n\n");
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ dihedral_class2_bb13_k[i],
+ dihedral_class2_bb13_r10[i],dihedral_class2_bb13_r30[i]);
+
+ } else if (strcmp(dihedral_style,"harmonic") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %d %d\n",i,
+ dihedral_harmonic_k[i],dihedral_harmonic_multiplicity[i],
+ dihedral_harmonic_sign[i]);
+
+ } else if (strcmp(dihedral_style,"cosineshiftexp") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,
+ dihedral_cosineshiftexp_umin[i],
+ dihedral_cosineshiftexp_theta[i]*180.0/PI,
+ dihedral_cosineshiftexp_a[i]
+ );
+
+ } else if (strcmp(dihedral_style,"helix") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,dihedral_helix_aphi[i],
+ dihedral_helix_bphi[i],dihedral_helix_cphi[i]);
+
+ } else if (strcmp(dihedral_style,"multi/harmonic") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g\n",i,
+ dihedral_multi_a1[i],dihedral_multi_a2[i],
+ dihedral_multi_a3[i],dihedral_multi_a4[i],
+ dihedral_multi_a5[i]);
+
+ } else if (strcmp(dihedral_style,"opls") == 0) {
+ for (int i = 1; i <= ndihedraltypes; i++) // restore factor of 2
+ fprintf(fp,"%d %g %g %g %g\n",i,
+ 2.0*dihedral_opls_k1[i],2.0*dihedral_opls_k2[i],
+ 2.0*dihedral_opls_k3[i],2.0*dihedral_opls_k4[i]);
+ }
+ }
+
+ if (improper_style) {
+ if ((strcmp(improper_style,"none") != 0) &&
+ (strcmp(improper_style,"hybrid") != 0))
+ fprintf(fp,"\nImproper Coeffs\n\n");
+
+ if (strcmp(improper_style,"class2") == 0) {
+ for (int i = 1; i <= nimpropertypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ improper_class2_k0[i],improper_class2_chi0[i]/PI*180.0);
+
+ fprintf(fp,"\nAngleAngle Coeffs\n\n");
+ for (int i = 1; i <= nimpropertypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g\n",i,
+ improper_class2_aa_k1[i],improper_class2_aa_k2[i],
+ improper_class2_aa_k3[i],
+ improper_class2_aa_theta0_1[i]/PI*180.0,
+ improper_class2_aa_theta0_2[i]/PI*180.0,
+ improper_class2_aa_theta0_3[i]/PI*180.0);
+
+ } else if (strcmp(improper_style,"cvff") == 0) {
+ for (int i = 1; i <= nimpropertypes; i++)
+ fprintf(fp,"%d %g %d %d\n",i,
+ improper_cvff_k[i],improper_cvff_sign[i],
+ improper_cvff_multiplicity[i]);
+
+ } else if (strcmp(improper_style,"harmonic") == 0) {
+ for (int i = 1; i <= nimpropertypes; i++)
+ fprintf(fp,"%d %g %g\n",i,
+ improper_harmonic_k[i],improper_harmonic_chi[i]/PI*180.0);
+ }
+ }
+
+ if (natoms) {
+ fprintf(fp,"\nAtoms\n\n");
+
+ int ix,iy,iz;
+ for (bigint i = 0; i < natoms; i++) {
+
+ ix = (image[i] & 1023) - 512;
+ iy = (image[i] >> 10 & 1023) - 512;
+ iz = (image[i] >> 20) - 512;
+
+ if (style_hybrid == 0) {
+ if (style_angle) write_atom_angle(fp,i,ix,iy,iz);
+ if (style_atomic) write_atom_atomic(fp,i,ix,iy,iz);
+ if (style_bond) write_atom_bond(fp,i,ix,iy,iz);
+ if (style_charge) write_atom_charge(fp,i,ix,iy,iz);
+ if (style_dipole) write_atom_dipole(fp,i,ix,iy,iz);
+ if (style_ellipsoid) write_atom_ellipsoid(fp,i,ix,iy,iz);
+ if (style_full) write_atom_full(fp,i,ix,iy,iz);
+ if (style_sphere) write_atom_sphere(fp,i,ix,iy,iz);
+ if (style_molecular) write_atom_molecular(fp,i,ix,iy,iz);
+ if (style_peri) write_atom_peri(fp,i,ix,iy,iz);
+ fprintf(fp,"\n");
+
+ } else {
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e",
+ tag[i],type[i],x[i],y[i],z[i]);
+ for (int k = 1; k <= style_hybrid; k++) {
+ if (k == style_angle) write_atom_angle_extra(fp,i);
+ if (k == style_atomic) write_atom_atomic_extra(fp,i);
+ if (k == style_bond) write_atom_bond_extra(fp,i);
+ if (k == style_charge) write_atom_charge_extra(fp,i);
+ if (k == style_dipole) write_atom_dipole_extra(fp,i);
+ if (k == style_ellipsoid) write_atom_ellipsoid_extra(fp,i);
+ if (k == style_full) write_atom_full_extra(fp,i);
+ if (k == style_sphere) write_atom_sphere_extra(fp,i);
+ if (k == style_molecular) write_atom_molecular_extra(fp,i);
+ if (k == style_peri) write_atom_peri_extra(fp,i);
+ }
+ fprintf(fp," %d %d %d\n",ix,iy,iz);
+ }
+ }
+ }
+
+ if (natoms) {
+ fprintf(fp,"\nVelocities\n\n");
+ for (bigint i = 0; i < natoms; i++)
+
+ if (style_hybrid == 0) {
+ if (style_angle) write_vel_angle(fp,i);
+ if (style_atomic) write_vel_atomic(fp,i);
+ if (style_bond) write_vel_bond(fp,i);
+ if (style_charge) write_vel_charge(fp,i);
+ if (style_dipole) write_vel_dipole(fp,i);
+ if (style_ellipsoid) write_vel_ellipsoid(fp,i);
+ if (style_full) write_vel_full(fp,i);
+ if (style_sphere) write_vel_sphere(fp,i);
+ if (style_molecular) write_vel_molecular(fp,i);
+ if (style_peri) write_vel_peri(fp,i);
+ fprintf(fp,"\n");
+
+ } else {
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+ for (int k = 1; k <= style_hybrid; k++) {
+ if (k == style_angle) write_vel_angle_extra(fp,i);
+ if (k == style_atomic) write_vel_atomic_extra(fp,i);
+ if (k == style_bond) write_vel_bond_extra(fp,i);
+ if (k == style_charge) write_vel_charge_extra(fp,i);
+ if (k == style_dipole) write_vel_dipole_extra(fp,i);
+ if (k == style_ellipsoid) write_vel_ellipsoid_extra(fp,i);
+ if (k == style_full) write_vel_full_extra(fp,i);
+ if (k == style_sphere) write_vel_sphere_extra(fp,i);
+ if (k == style_molecular) write_vel_molecular_extra(fp,i);
+ if (k == style_peri) write_vel_peri_extra(fp,i);
+ }
+ fprintf(fp,"\n");
+ }
+ }
+
+ if (nellipsoids) {
+ fprintf(fp,"\nEllipsoids\n\n");
+ for (bigint i = 0; i < natoms; i++) {
+ if (ellipsoid[i])
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e "
+ "%-1.16e %-1.16e %-1.16e %-1.16e \n",
+ tag[i],2.0*shapex[i],2.0*shapey[i],2.0*shapez[i],
+ quatw[i],quati[i],quatj[i],quatk[i]);
+ }
+ }
+
+ if (nbonds) {
+ fprintf(fp,"\nBonds\n\n");
+ for (bigint i = 0; i < nbonds; i++)
+ fprintf(fp,BIGINT_FORMAT " %d %d %d\n",
+ i+1,bond_type[i],bond_atom1[i],bond_atom2[i]);
+ }
+
+ if (nangles) {
+ fprintf(fp,"\nAngles\n\n");
+ for (bigint i = 0; i < nangles; i++)
+ fprintf(fp,BIGINT_FORMAT " %d %d %d %d\n",
+ i+1,angle_type[i],angle_atom1[i],angle_atom2[i],angle_atom3[i]);
+ }
+
+ if (ndihedrals) {
+ fprintf(fp,"\nDihedrals\n\n");
+ for (bigint i = 0; i < ndihedrals; i++)
+ fprintf(fp,BIGINT_FORMAT " %d %d %d %d %d\n",
+ i+1,dihedral_type[i],dihedral_atom1[i],dihedral_atom2[i],
+ dihedral_atom3[i],dihedral_atom4[i]);
+ }
+
+ if (nimpropers) {
+ fprintf(fp,"\nImpropers\n\n");
+ for (bigint i = 0; i < nimpropers; i++)
+ fprintf(fp,BIGINT_FORMAT " %d %d %d %d %d\n",
+ i+1,improper_type[i],improper_atom1[i],improper_atom2[i],
+ improper_atom3[i],improper_atom4[i]);
+ }
+}
+
+// ---------------------------------------------------------------------
+// per-atom write routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_atom_angle(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_atomic(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_bond(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_charge(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],q[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_dipole(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
+ "%-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],q[i],x[i],y[i],z[i],
+ mux[i],muy[i],muz[i],ix,iy,iz);
+}
+
+void Data::write_atom_ellipsoid(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],ellipsoid[i],density[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_full(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],molecule[i],type[i],q[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_sphere(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],2.0*radius[i],density[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_molecular(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_peri(FILE *fp, int i, int ix, int iy, int iz)
+{
+ fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+ tag[i],type[i],vfrac[i],rmass[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+// ---------------------------------------------------------------------
+// per-atom write routines of extra quantities unique to style
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_atom_angle_extra(FILE *fp, int i)
+{
+ fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_atomic_extra(FILE *fp, int i) {}
+
+void Data::write_atom_bond_extra(FILE *fp, int i)
+{
+ fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_charge_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e",q[i]);
+}
+
+void Data::write_atom_dipole_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",q[i],mux[i],muy[i],muz[i]);
+}
+
+void Data::write_atom_ellipsoid_extra(FILE *fp, int i)
+{
+ fprintf(fp," %d %-1.16e",ellipsoid[i],density[i]);
+}
+
+void Data::write_atom_full_extra(FILE *fp, int i)
+{
+ fprintf(fp," %d %-1.16e",molecule[i],q[i]);
+}
+
+void Data::write_atom_sphere_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e %-1.16e",2.0*radius[i],density[i]);
+}
+
+void Data::write_atom_molecular_extra(FILE *fp, int i)
+{
+ fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_peri_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e %-1.16e",vfrac[i],rmass[i]);
+}
+
+// ---------------------------------------------------------------------
+// per-atom velocity write routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_vel_angle(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_atomic(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_bond(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_charge(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_dipole(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_ellipsoid(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
+ tag[i],vx[i],vy[i],vz[i],angmomx[i],angmomy[i],angmomz[i]);
+}
+
+void Data::write_vel_full(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_sphere(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
+ tag[i],vx[i],vy[i],vz[i],omegax[i],omegay[i],omegaz[i]);
+}
+
+void Data::write_vel_molecular(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_peri(FILE *fp, int i)
+{
+ fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+// ---------------------------------------------------------------------
+// per-atom velocity write routines of extra quantities unique to style
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_vel_angle_extra(FILE *fp, int i) {}
+void Data::write_vel_atomic_extra(FILE *fp, int i) {}
+void Data::write_vel_bond_extra(FILE *fp, int i) {}
+void Data::write_vel_charge_extra(FILE *fp, int i) {}
+void Data::write_vel_dipole_extra(FILE *fp, int i) {}
+
+void Data::write_vel_ellipsoid_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e %-1.16e %-1.16e",angmomx[i],angmomy[i],angmomz[i]);
+}
+
+void Data::write_vel_full_extra(FILE *fp, int i) {}
+
+void Data::write_vel_sphere_extra(FILE *fp, int i)
+{
+ fprintf(fp," %-1.16e %-1.16e %-1.16e",omegax[i],omegay[i],omegaz[i]);
+}
+
+void Data::write_vel_molecular_extra(FILE *fp, int i) {}
+void Data::write_vel_peri_extra(FILE *fp, int i) {}
+
+// ---------------------------------------------------------------------
+// strip known accelerator suffixes from style name
+// ---------------------------------------------------------------------
+
+void strip_suffix(char *style)
+{
+ char *slash = strrchr(style,'/');
+ if (slash == NULL) return;
+
+ int i=0;
+
+ const char *suffix_list[] = { "/opt", "/gpu", "/cuda", "/omp", NULL };
+ const char *suffix = suffix_list[0];
+ while (suffix != NULL) {
+ if (strcmp(slash,suffix) == 0) {
+ *slash = '\0';
+ return;
+ }
+ ++i;
+ suffix = suffix_list[i];
+ }
+}
+
+// ---------------------------------------------------------------------
+// binary reads from restart file
+// ---------------------------------------------------------------------
+
+int read_int(FILE *fp)
+{
+ int value;
+ fread(&value,sizeof(int),1,fp);
+ return value;
+}
+
+double read_double(FILE *fp)
+{
+ double value;
+ fread(&value,sizeof(double),1,fp);
+ return value;
+}
+
+char *read_char(FILE *fp)
+{
+ int n;
+ fread(&n,sizeof(int),1,fp);
+ if (n == 0) return NULL;
+ char *value = new char[n];
+ fread(value,sizeof(char),n,fp);
+ return value;
+}
+
+bigint read_bigint(FILE *fp)
+{
+ bigint value;
+ fread(&value,sizeof(bigint),1,fp);
+ return value;
+}
More information about the Swift-commit
mailing list