Bug fix in OpenMM file generation

.gitignore

@@ -1,8 +1,3 @@
-examples/*/*.xyz
-examples/*/*.lmp
-examples/*/*.lammpstrj
-examples/*/*.log
-examples/*/*.pdb
 bin/
 docs/html/
 lib/libpackmol.a

docs/Doxyfile

@@ -1,6 +1,6 @@
 DOXYFILE_ENCODING      = UTF-8
 PROJECT_NAME           = "Playmol"
-PROJECT_BRIEF = "(Version: 25 Mar 2018)"
+PROJECT_BRIEF = "(Version: 11 Apr 2018)"
 PROJECT_LOGO           = playmol_logo.gif
 OUTPUT_DIRECTORY       = .
 OUTPUT_DIRECTORY       = 

examples/.gitignore

@@ -0,0 +1,7 @@
+log.lammps
+playmol.log
+*.xml
+*.pdb
+*.lammpstrj
+*.lmp
+*.xyz

examples/tip4p/openmm.py

@@ -6,7 +6,7 @@
 
 model = 'tip4p-2005'
 resorces = {1:'CPU', 2:'CUDA', 3:'OpenCL'}
-choice = 1
+choice = 2
 temp = 298*unit.kelvin
 press = 1.0*unit.atmosphere
 dt = 1.0*unit.femtosecond
@@ -42,4 +42,5 @@
 simulation.context.setVelocitiesToTemperature(temp)
 
 print('Running Production...')
-simulation.step(nsteps)
+for i in range(nsteps):
+    simulation.step(1)

examples/tip4p/tip4p.playmol

@@ -21,9 +21,11 @@ diameter     OW 3.1
 diameter     HW 2.0
 diameter     MW 0.1
 
-bond_type    OW HW zero $rOH
+bond_type    OW HW harmonic 553 $rOH
 bond_type    OW MW zero $rOM
 
+angle_type   HW OW HW harmonic 100 $theta0
+
 atom         OW  OW  0.0000
 atom         HW1 HW  $qH
 atom         HW2 HW  $qH

src/playmol.f90

@@ -24,7 +24,7 @@ program playmol
 integer :: i, inp
 character(sl) :: infile
 type(tPlaymol) :: System
-call writeln( "Playmol (Version: 25 Mar 2018)" )
+call writeln( "Playmol (Version: 11 Apr 2018)" )
 if (iargc() == 0) call error( "Usage: playmol <file-1> <file-2> ..." )
 call init_log( file = "playmol.log" )
 call getarg( 1, infile )

src/write_lammps.f90

@@ -254,13 +254,17 @@ end subroutine write_type
       !---------------------------------------------------------------------------------------------
       subroutine write_masses
         integer :: i
-        character(sl) :: mass, element
+        character(sl) :: mass, E(size(atom_types))
         write(unit,'(/,"Masses",/)')
         do i = 1, size(atom_types)
-          call me % element_and_mass( atom_types(i)%types, element, mass )
+          call me % element_and_mass( atom_types(i)%types, E(i), mass )
           if (mass == real2str(0.0_rb)) mass = "1.0E-20"
           write(unit,'(A)') trim(join([int2str(i), mass, "#", atom_types(i)%types]))
+          if (E(i) == "UA") then
+            E(i) = trim(me%elements(minloc(abs(me%masses - str2real(mass)), dim = 1)))//"?"
+          end if
         end do
+        write(unit,'(/,"# Elements: ",A)') trim(join(E))
       end subroutine write_masses
       !---------------------------------------------------------------------------------------------
       subroutine write_atoms( natoms )

src/write_openmm.f90

@@ -68,6 +68,7 @@ subroutine tPlaymol_write_openmm( me, unit, keywords )
       end do
     end block
 
+    write(unit,'("<!--Generated by Playmol on ",A,"-->")') trim(now())
     write(unit,'("<ForceField>")')
     call atom_types()
 
@@ -289,7 +290,7 @@ subroutine bond_types()
             call error( "harmonic bond model required" )
           end if
           call items(4, "Bond", p, [current%id, real2str([r0, K])])
-        current => current % next
+          current => current % next
         end do
         call list % destroy(silent = .true.)
         write(unit,'(2X,"</HarmonicBondForce>")')
@@ -383,11 +384,12 @@ subroutine non_bonded_model()
         type(StrucList) :: list
         type(Struc), pointer :: current
         real(rb), parameter :: tol = 1.0E-8_rb
-        character(sl), parameter :: zero = "0", &
-                                    p(4) = [character(sl) :: "type", "charge", "sigma", "epsilon"]
+        character(sl), parameter :: zero = "0", one = "1", &
+                                    p(3) = [character(sl) :: "type", "sigma", "epsilon"]
 
         write(unit,'(2X,"<NonbondedForce ",A,X,A,">")') trim(item("coulomb14scale", coul14)), &
                                                         trim(item("lj14scale", lj14))
+        write(unit,'(4X,"<UseAttributeFromResidue name=""charge""/>")')
         list = local_list( me % atom_list, me % atom_type_list, .false. )
         current => list % first
         do while (associated(current))
@@ -400,11 +402,11 @@ subroutine non_bonded_model()
             sig = real2str(str2real(arg(3)) * length)
           else if ((arg(1) == "zero").or.(arg(1)(1:4) == "coul")) then
             eps = zero
-            sig = zero
+            sig = one
           else
             call error( "Lennard-Jones potential model required" )
           end if
-          call items(4, "Atom", p, [current%id(1), zero, sig, eps])
+          call items(4, "Atom", p, [current%id(1), sig, eps])
           current => current % next
         end do
         call list % destroy(silent = .true.)

src/write_pdb.f90

@@ -50,7 +50,6 @@ subroutine tPlaymol_write_pdb( me, unit )
     do while (associated(current))
       jmol = jmol + 1
       imol = str2int(me % molecules % list % parameters( current % id ) )
-!      residue = me % molecule_names( imol )merge(HOH, letterCode(imol - count(water(1:imol-1))), water(imol))
       do i = 1, natoms(imol)
         iatom = iatom + 1
         call me % atom_list % search( current%id, ptr, indx )
@@ -69,7 +68,7 @@ subroutine tPlaymol_write_pdb( me, unit )
           element ! Element symbol
         current => current % next
       end do
-!      if (.not.water(imol)) write(unit,'("TER")')
+      write(unit,'("TER")')
     end do
 
     ! Write out connectivity info: