* generation of previously equilibrated SDS micelle
*

!identifiers
set id1 micelle
set id2 1

! 1: topology
! 2: parameter
! 3: cutim
! 4: cutnb
! 5: ctonnb
! 6: ctofnb
! 7: electrostatic cutoff regime (shift or switch)
! 8: electrostatic pairing method (atom or group)
! 9: VDW pairing method (vatom or vgroup)
! VDW cutoff regime is vwsitch
set 1 top_all27_lipid.rtf
set 2 par_all27_lipid.prm
set 3 12.0 ! cutim
set 4 12.0 ! cutnb
set 5  8.0 ! ctonnb
set 6 10.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

! read topology and parameters
open unit 9 read form name @1
read rtf card unit 9

open unit 9 read form name @2
read para card unit 9

!generate SDS
read sequence SDS 60

generate mice first none last none setup warn

!generate water
read sequence TIP3 4458

generate wat first none last none setup warn noangle nodihedral

!generate sodium
read sequence SOD 60

generate sod first none last none setup warn noangle nodihedral

!read coordinates
open unit 20 read form name dyn_mic4_box_120ps.crd
read coor card unit 20 resid

!delete uninitialized (no coordinates) waters; these correspond
!to the 60 waters that were initially deleted and replaced
!with sodiums
delete atom sele .not. init show end

!write psf and coordinates
open unit 20 write form name @id1_@id2.psf
write psf card unit 20
* Previously equilibrated 60 SDS Micelle with water and sodium
* XDIM = 54.14101  YDIM = 54.14101  ZDIM = 54.14101
*

open unit 20 write form name @id1_@id2.crd
write coor card unit 20
* Previously equilibrated 60 SDS Micelle with water and sodium
* XDIM = 54.14101  YDIM = 54.14101  ZDIM = 54.14101
*

! cell parameters
set a 54.14101
set b 54.14101
set c 54.14101 

crystal define cubic @a @b @c 90. 90. 90.

crystal build noper 0 cutoff @3

!turn on image centering of solvent
image byres xcen 0.0 ycen 0.0 zcen 0.0 -
     sele resn tip3 .or. resn sod end

!kappa should be optimized based on the value of
!kappa at which the 1st derivative of the energy
!with respect to kappa goes to 0.0
set kap 0.36  !NOT Optimized
set ord 6

!fftx dimensions should correspond to 1 A or less; must
!be power of 2, 3 or 5
energy inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx 64 ffty 64 fftz 64 order @ord -
@7 @8 @9 vfswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5

stop