8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
874 if (oldion.eq.1) then
876 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
877 print *,msc(i,5),restok(i,5)
881 read (iion,*) ncatprotparm
882 allocate(catprm(ncatprotparm,4))
884 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
888 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
892 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
895 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
896 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
900 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
901 restyp(i,5),"-",restyp(j,2)
904 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
905 !----------------------------------------------------
906 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
907 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
908 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
909 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
910 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
911 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
921 allocate(liptranene(ntyp))
922 !C reading lipid parameters
923 write (iout,*) "iliptranpar",iliptranpar
925 read(iliptranpar,*) pepliptran
928 read(iliptranpar,*) liptranene(i)
929 print *,liptranene(i)
935 ! Read the parameters of the probability distribution/energy expression
936 ! of the virtual-bond valence angles theta
939 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
940 (bthet(j,i,1,1),j=1,2)
941 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
942 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
943 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
947 athet(1,i,1,-1)=athet(1,i,1,1)
948 athet(2,i,1,-1)=athet(2,i,1,1)
949 bthet(1,i,1,-1)=-bthet(1,i,1,1)
950 bthet(2,i,1,-1)=-bthet(2,i,1,1)
951 athet(1,i,-1,1)=-athet(1,i,1,1)
952 athet(2,i,-1,1)=-athet(2,i,1,1)
953 bthet(1,i,-1,1)=bthet(1,i,1,1)
954 bthet(2,i,-1,1)=bthet(2,i,1,1)
958 athet(1,i,-1,-1)=athet(1,-i,1,1)
959 athet(2,i,-1,-1)=-athet(2,-i,1,1)
960 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
961 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
962 athet(1,i,-1,1)=athet(1,-i,1,1)
963 athet(2,i,-1,1)=-athet(2,-i,1,1)
964 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
965 bthet(2,i,-1,1)=bthet(2,-i,1,1)
966 athet(1,i,1,-1)=-athet(1,-i,1,1)
967 athet(2,i,1,-1)=athet(2,-i,1,1)
968 bthet(1,i,1,-1)=bthet(1,-i,1,1)
969 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
974 polthet(j,i)=polthet(j,-i)
977 gthet(j,i)=gthet(j,-i)
985 'Parameters of the virtual-bond valence angles:'
986 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
987 ' ATHETA0 ',' A1 ',' A2 ',&
990 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
991 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
993 write (iout,'(/a/9x,5a/79(1h-))') &
994 'Parameters of the expression for sigma(theta_c):',&
995 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
996 ' ALPH3 ',' SIGMA0C '
998 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
999 (polthet(j,i),j=0,3),sigc0(i)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the second gaussian:',&
1003 ' THETA0 ',' SIGMA0 ',' G1 ',&
1006 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1007 sig0(i),(gthet(j,i),j=1,3)
1010 write (iout,'(a)') &
1011 'Parameters of the virtual-bond valence angles:'
1012 write (iout,'(/a/9x,5a/79(1h-))') &
1013 'Coefficients of expansion',&
1014 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1015 ' b1*10^1 ',' b2*10^1 '
1017 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1018 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1019 (10*bthet(j,i,1,1),j=1,2)
1021 write (iout,'(/a/9x,5a/79(1h-))') &
1022 'Parameters of the expression for sigma(theta_c):',&
1023 ' alpha0 ',' alph1 ',' alph2 ',&
1024 ' alhp3 ',' sigma0c '
1026 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1027 (polthet(j,i),j=0,3),sigc0(i)
1029 write (iout,'(/a/9x,5a/79(1h-))') &
1030 'Parameters of the second gaussian:',&
1031 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1034 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1035 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1041 ! Read the parameters of Utheta determined from ab initio surfaces
1042 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1044 IF (tor_mode.eq.0) THEN
1045 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1046 ntheterm3,nsingle,ndouble
1047 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1049 !----------------------------------------------------
1050 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1051 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1052 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1053 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1054 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1055 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1056 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1057 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1058 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1059 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1060 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1061 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1062 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1063 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1064 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1065 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1066 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1067 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1068 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1069 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1070 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1071 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1072 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1073 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1075 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1077 ithetyp(i)=-ithetyp(-i)
1080 aa0thet(:,:,:,:)=0.0d0
1081 aathet(:,:,:,:,:)=0.0d0
1082 bbthet(:,:,:,:,:,:)=0.0d0
1083 ccthet(:,:,:,:,:,:)=0.0d0
1084 ddthet(:,:,:,:,:,:)=0.0d0
1085 eethet(:,:,:,:,:,:)=0.0d0
1086 ffthet(:,:,:,:,:,:,:)=0.0d0
1087 ggthet(:,:,:,:,:,:,:)=0.0d0
1089 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1091 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1092 ! VAR:1=non-glicyne non-proline 2=proline
1093 ! VAR:negative values for D-aminoacid
1095 do j=-nthetyp,nthetyp
1096 do k=-nthetyp,nthetyp
1097 read (ithep,'(6a)',end=111,err=111) res1
1098 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1099 ! VAR: aa0thet is variable describing the average value of Foureir
1100 ! VAR: expansion series
1101 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1102 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1103 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1104 read (ithep,*,end=111,err=111) &
1105 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1106 read (ithep,*,end=111,err=111) &
1107 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1108 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1109 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1110 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1112 read (ithep,*,end=111,err=111) &
1113 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1114 ffthet(lll,llll,ll,i,j,k,iblock),&
1115 ggthet(llll,lll,ll,i,j,k,iblock),&
1116 ggthet(lll,llll,ll,i,j,k,iblock),&
1117 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1122 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1123 ! coefficients of theta-and-gamma-dependent terms are zero.
1124 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1125 ! RECOMENTDED AFTER VERSION 3.3)
1129 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1130 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1132 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1133 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1136 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1138 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1141 ! AND COMMENT THE LOOPS BELOW
1145 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1146 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1148 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1149 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1152 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1154 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1159 ! Substitution for D aminoacids from symmetry.
1162 do j=-nthetyp,nthetyp
1163 do k=-nthetyp,nthetyp
1164 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1166 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1170 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1171 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1172 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1173 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1179 ffthet(llll,lll,ll,i,j,k,iblock)= &
1180 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1181 ffthet(lll,llll,ll,i,j,k,iblock)= &
1182 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1183 ggthet(llll,lll,ll,i,j,k,iblock)= &
1184 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1185 ggthet(lll,llll,ll,i,j,k,iblock)= &
1186 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1195 ! Control printout of the coefficients of virtual-bond-angle potentials
1198 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1203 write (iout,'(//4a)') &
1204 'Type ',onelett(i),onelett(j),onelett(k)
1205 write (iout,'(//a,10x,a)') " l","a[l]"
1206 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1207 write (iout,'(i2,1pe15.5)') &
1208 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1210 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1211 "b",l,"c",l,"d",l,"e",l
1213 write (iout,'(i2,4(1pe15.5))') m,&
1214 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1215 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1219 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1220 "f+",l,"f-",l,"g+",l,"g-",l
1223 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1224 ffthet(n,m,l,i,j,k,iblock),&
1225 ffthet(m,n,l,i,j,k,iblock),&
1226 ggthet(n,m,l,i,j,k,iblock),&
1227 ggthet(m,n,l,i,j,k,iblock)
1238 !C here will be the apropriate recalibrating for D-aminoacid
1239 read (ithep,*,end=121,err=121) nthetyp
1240 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1241 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1242 do i=-nthetyp+1,nthetyp-1
1243 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1244 do j=0,nbend_kcc_Tb(i)
1245 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1249 write (iout,'(a)') &
1250 "Parameters of the valence-only potentials"
1251 do i=-nthetyp+1,nthetyp-1
1252 write (iout,'(2a)') "Type ",toronelet(i)
1253 do k=0,nbend_kcc_Tb(i)
1254 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1260 write (2,*) "Start reading THETA_PDB",ithep_pdb
1262 ! write (2,*) 'i=',i
1263 read (ithep_pdb,*,err=111,end=111) &
1264 a0thet(i),(athet(j,i,1,1),j=1,2),&
1265 (bthet(j,i,1,1),j=1,2)
1266 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1267 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1268 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1269 sigc0(i)=sigc0(i)**2
1272 athet(1,i,1,-1)=athet(1,i,1,1)
1273 athet(2,i,1,-1)=athet(2,i,1,1)
1274 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1275 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1276 athet(1,i,-1,1)=-athet(1,i,1,1)
1277 athet(2,i,-1,1)=-athet(2,i,1,1)
1278 bthet(1,i,-1,1)=bthet(1,i,1,1)
1279 bthet(2,i,-1,1)=bthet(2,i,1,1)
1282 a0thet(i)=a0thet(-i)
1283 athet(1,i,-1,-1)=athet(1,-i,1,1)
1284 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1285 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1286 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1287 athet(1,i,-1,1)=athet(1,-i,1,1)
1288 athet(2,i,-1,1)=-athet(2,-i,1,1)
1289 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1290 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1291 athet(1,i,1,-1)=-athet(1,-i,1,1)
1292 athet(2,i,1,-1)=athet(2,-i,1,1)
1293 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1294 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1295 theta0(i)=theta0(-i)
1299 polthet(j,i)=polthet(j,-i)
1302 gthet(j,i)=gthet(j,-i)
1305 write (2,*) "End reading THETA_PDB"
1309 !--------------- Reading theta parameters for nucleic acid-------
1310 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1311 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1312 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1313 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1314 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1315 nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1317 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1318 nthetyp_nucl+1,nthetyp_nucl+1))
1319 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1320 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1321 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1322 nthetyp_nucl+1,nthetyp_nucl+1))
1323 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1324 nthetyp_nucl+1,nthetyp_nucl+1))
1325 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1326 nthetyp_nucl+1,nthetyp_nucl+1))
1327 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1328 nthetyp_nucl+1,nthetyp_nucl+1))
1329 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1330 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1331 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1332 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1333 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1334 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1336 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1337 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1339 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1341 aa0thet_nucl(:,:,:)=0.0d0
1342 aathet_nucl(:,:,:,:)=0.0d0
1343 bbthet_nucl(:,:,:,:,:)=0.0d0
1344 ccthet_nucl(:,:,:,:,:)=0.0d0
1345 ddthet_nucl(:,:,:,:,:)=0.0d0
1346 eethet_nucl(:,:,:,:,:)=0.0d0
1347 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1348 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1353 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1354 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1355 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1356 read (ithep_nucl,*,end=111,err=111) &
1357 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1358 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1359 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1360 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1361 read (ithep_nucl,*,end=111,err=111) &
1362 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1363 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1364 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1369 !-------------------------------------------
1370 allocate(nlob(ntyp1)) !(ntyp1)
1371 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1372 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1373 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1380 gaussc(:,:,:,:)=0.0D0
1384 ! Read the parameters of the probability distribution/energy expression
1385 ! of the side chains.
1388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1392 dsc_inv(i)=1.0D0/dsc(i)
1403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1404 ((blower(k,l,1),l=1,k),k=1,3)
1405 censc(1,1,-i)=censc(1,1,i)
1406 censc(2,1,-i)=censc(2,1,i)
1407 censc(3,1,-i)=-censc(3,1,i)
1409 read (irotam,*,end=112,err=112) bsc(j,i)
1410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1411 ((blower(k,l,j),l=1,k),k=1,3)
1412 censc(1,j,-i)=censc(1,j,i)
1413 censc(2,j,-i)=censc(2,j,i)
1414 censc(3,j,-i)=-censc(3,j,i)
1415 ! BSC is amplitude of Gaussian
1422 akl=akl+blower(k,m,j)*blower(l,m,j)
1426 if (((k.eq.3).and.(l.ne.3)) &
1427 .or.((l.eq.3).and.(k.ne.3))) then
1428 gaussc(k,l,j,-i)=-akl
1429 gaussc(l,k,j,-i)=-akl
1431 gaussc(k,l,j,-i)=akl
1432 gaussc(l,k,j,-i)=akl
1441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1444 if (nlobi.gt.0) then
1446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1447 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1449 'log h',(bsc(j,i),j=1,nlobi)
1450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1451 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1454 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1458 write (iout,'(a,f10.4,4(16x,f10.4))') &
1459 'Center ',(bsc(j,i),j=1,nlobi)
1460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1469 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1470 ! added by Urszula Kozlowska 07/11/2007
1472 !el Maximum number of SC local term fitting function coefficiants
1473 !el integer,parameter :: maxsccoef=65
1475 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1478 read (irotam,*,end=112,err=112)
1480 read (irotam,*,end=112,err=112)
1483 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1487 !---------reading nucleic acid parameters for rotamers-------------------
1488 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1489 do i=1,ntyp_molec(2)
1490 read (irotam_nucl,*,end=112,err=112)
1492 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1498 write (iout,*) "Base rotamer parameters"
1499 do i=1,ntyp_molec(2)
1500 write (iout,'(a)') restyp(i,2)
1501 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1506 ! Read the parameters of the probability distribution/energy expression
1507 ! of the side chains.
1509 write (2,*) "Start reading ROTAM_PDB"
1511 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1515 dsc_inv(i)=1.0D0/dsc(i)
1526 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1527 ((blower(k,l,1),l=1,k),k=1,3)
1529 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1530 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1531 ((blower(k,l,j),l=1,k),k=1,3)
1538 akl=akl+blower(k,m,j)*blower(l,m,j)
1548 write (2,*) "End reading ROTAM_PDB"
1554 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1555 !C interaction energy of the Gly, Ala, and Pro prototypes.
1557 read (ifourier,*) nloctyp
1558 SPLIT_FOURIERTOR = nloctyp.lt.0
1559 nloctyp = iabs(nloctyp)
1560 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1561 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1562 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1563 !C allocate(ctilde(2,2,nres))
1564 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1565 !C allocate(gtb1(2,nres))
1566 !C allocate(gtb2(2,nres))
1567 !C allocate(cc(2,2,nres))
1568 !C allocate(dd(2,2,nres))
1569 !C allocate(ee(2,2,nres))
1570 !C allocate(gtcc(2,2,nres))
1571 !C allocate(gtdd(2,2,nres))
1572 !C allocate(gtee(2,2,nres))
1575 allocate(itype2loc(-ntyp1:ntyp1))
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(bnew1(3,2,-nloctyp:nloctyp))
1578 allocate(bnew2(3,2,-nloctyp:nloctyp))
1579 allocate(ccnew(3,2,-nloctyp:nloctyp))
1580 allocate(ddnew(3,2,-nloctyp:nloctyp))
1581 allocate(e0new(3,-nloctyp:nloctyp))
1582 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1583 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1584 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1585 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1586 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1587 allocate(e0newtor(3,-nloctyp:nloctyp))
1588 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1590 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1591 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1592 itype2loc(ntyp1)=nloctyp
1593 iloctyp(nloctyp)=ntyp1
1595 itype2loc(-i)=-itype2loc(i)
1598 allocate(iloctyp(-nloctyp:nloctyp))
1599 allocate(itype2loc(-ntyp1:ntyp1))
1606 iloctyp(-i)=-iloctyp(i)
1608 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1609 !c write (iout,*) "nloctyp",nloctyp,
1610 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1611 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1612 !c write (iout,*) "nloctyp",nloctyp,
1613 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1616 !c write (iout,*) "NEWCORR",i
1617 read (ifourier,*,end=115,err=115)
1620 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1623 !c write (iout,*) "NEWCORR BNEW1"
1624 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1630 !c write (iout,*) "NEWCORR BNEW2"
1631 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1633 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1634 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1636 !c write (iout,*) "NEWCORR CCNEW"
1637 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1639 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1640 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1642 !c write (iout,*) "NEWCORR DDNEW"
1643 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1647 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1651 !c write (iout,*) "NEWCORR EENEW1"
1652 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1654 read (ifourier,*,end=115,err=115) e0new(ii,i)
1656 !c write (iout,*) (e0new(ii,i),ii=1,3)
1658 !c write (iout,*) "NEWCORR EENEW"
1661 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1662 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1663 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1664 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1669 bnew1(ii,1,-i)= bnew1(ii,1,i)
1670 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1671 bnew2(ii,1,-i)= bnew2(ii,1,i)
1672 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1675 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1676 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1677 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1678 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1679 ccnew(ii,1,-i)=ccnew(ii,1,i)
1680 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1681 ddnew(ii,1,-i)=ddnew(ii,1,i)
1682 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1684 e0new(1,-i)= e0new(1,i)
1685 e0new(2,-i)=-e0new(2,i)
1686 e0new(3,-i)=-e0new(3,i)
1688 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1689 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1690 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1691 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1695 write (iout,'(a)') "Coefficients of the multibody terms"
1696 do i=-nloctyp+1,nloctyp-1
1697 write (iout,*) "Type: ",onelet(iloctyp(i))
1698 write (iout,*) "Coefficients of the expansion of B1"
1700 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1702 write (iout,*) "Coefficients of the expansion of B2"
1704 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1706 write (iout,*) "Coefficients of the expansion of C"
1707 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1708 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1709 write (iout,*) "Coefficients of the expansion of D"
1710 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1711 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1712 write (iout,*) "Coefficients of the expansion of E"
1713 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1716 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1721 IF (SPLIT_FOURIERTOR) THEN
1723 !c write (iout,*) "NEWCORR TOR",i
1724 read (ifourier,*,end=115,err=115)
1727 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1730 !c write (iout,*) "NEWCORR BNEW1 TOR"
1731 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1737 !c write (iout,*) "NEWCORR BNEW2 TOR"
1738 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1740 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1741 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1743 !c write (iout,*) "NEWCORR CCNEW TOR"
1744 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1746 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1747 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1749 !c write (iout,*) "NEWCORR DDNEW TOR"
1750 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1754 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1758 !c write (iout,*) "NEWCORR EENEW1 TOR"
1759 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1761 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1763 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1765 !c write (iout,*) "NEWCORR EENEW TOR"
1768 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1769 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1770 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1771 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1776 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1777 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1778 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1779 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1782 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1783 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1784 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1785 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1787 e0newtor(1,-i)= e0newtor(1,i)
1788 e0newtor(2,-i)=-e0newtor(2,i)
1789 e0newtor(3,-i)=-e0newtor(3,i)
1791 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1792 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1793 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1794 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1798 write (iout,'(a)') &
1799 "Single-body coefficients of the torsional potentials"
1800 do i=-nloctyp+1,nloctyp-1
1801 write (iout,*) "Type: ",onelet(iloctyp(i))
1802 write (iout,*) "Coefficients of the expansion of B1tor"
1804 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1805 j,(bnew1tor(k,j,i),k=1,3)
1807 write (iout,*) "Coefficients of the expansion of B2tor"
1809 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1810 j,(bnew2tor(k,j,i),k=1,3)
1812 write (iout,*) "Coefficients of the expansion of Ctor"
1813 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1814 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1815 write (iout,*) "Coefficients of the expansion of Dtor"
1816 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1817 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1818 write (iout,*) "Coefficients of the expansion of Etor"
1819 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1822 write (iout,'(1hE,2i1,2f10.5)') &
1823 j,k,(eenewtor(l,j,k,i),l=1,2)
1829 do i=-nloctyp+1,nloctyp-1
1832 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1837 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1841 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1842 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1843 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1844 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1847 ENDIF !SPLIT_FOURIER_TOR
1849 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1850 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1851 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1852 allocate(b(13,-nloctyp-1:nloctyp+1))
1854 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1856 read (ifourier,*,end=115,err=115)
1857 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1859 write (iout,*) 'Type ',onelet(iloctyp(i))
1860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1874 !c B1tilde(1,i) = b(3)
1875 !c! B1tilde(2,i) =-b(5)
1876 !c! B1tilde(1,-i) =-b(3)
1877 !c! B1tilde(2,-i) =b(5)
1878 !c! b1tilde(1,i)=0.0d0
1879 !c b1tilde(2,i)=0.0d0
1884 !cc B1tilde(1,i) = b(3,i)
1885 !cc B1tilde(2,i) =-b(5,i)
1886 !c B1tilde(1,-i) =-b(3,i)
1887 !c B1tilde(2,-i) =b(5,i)
1888 !cc b1tilde(1,i)=0.0d0
1889 !cc b1tilde(2,i)=0.0d0
1890 !cc B2(1,i) = b(2,i)
1891 !cc B2(2,i) = b(4,i)
1893 !c B2(2,-i) =-b(4,i)
1897 CCold(1,1,i)= b(7,i)
1898 CCold(2,2,i)=-b(7,i)
1899 CCold(2,1,i)= b(9,i)
1900 CCold(1,2,i)= b(9,i)
1901 CCold(1,1,-i)= b(7,i)
1902 CCold(2,2,-i)=-b(7,i)
1903 CCold(2,1,-i)=-b(9,i)
1904 CCold(1,2,-i)=-b(9,i)
1909 !c Ctilde(1,1,i)= CCold(1,1,i)
1910 !c Ctilde(1,2,i)= CCold(1,2,i)
1911 !c Ctilde(2,1,i)=-CCold(2,1,i)
1912 !c Ctilde(2,2,i)=-CCold(2,2,i)
1917 !c Ctilde(1,1,i)= CCold(1,1,i)
1918 !c Ctilde(1,2,i)= CCold(1,2,i)
1919 !c Ctilde(2,1,i)=-CCold(2,1,i)
1920 !c Ctilde(2,2,i)=-CCold(2,2,i)
1922 !c Ctilde(1,1,i)=0.0d0
1923 !c Ctilde(1,2,i)=0.0d0
1924 !c Ctilde(2,1,i)=0.0d0
1925 !c Ctilde(2,2,i)=0.0d0
1926 DDold(1,1,i)= b(6,i)
1927 DDold(2,2,i)=-b(6,i)
1928 DDold(2,1,i)= b(8,i)
1929 DDold(1,2,i)= b(8,i)
1930 DDold(1,1,-i)= b(6,i)
1931 DDold(2,2,-i)=-b(6,i)
1932 DDold(2,1,-i)=-b(8,i)
1933 DDold(1,2,-i)=-b(8,i)
1938 !c Dtilde(1,1,i)= DD(1,1,i)
1939 !c Dtilde(1,2,i)= DD(1,2,i)
1940 !c Dtilde(2,1,i)=-DD(2,1,i)
1941 !c Dtilde(2,2,i)=-DD(2,2,i)
1943 !c Dtilde(1,1,i)=0.0d0
1944 !c Dtilde(1,2,i)=0.0d0
1945 !c Dtilde(2,1,i)=0.0d0
1946 !c Dtilde(2,2,i)=0.0d0
1947 EEold(1,1,i)= b(10,i)+b(11,i)
1948 EEold(2,2,i)=-b(10,i)+b(11,i)
1949 EEold(2,1,i)= b(12,i)-b(13,i)
1950 EEold(1,2,i)= b(12,i)+b(13,i)
1951 EEold(1,1,-i)= b(10,i)+b(11,i)
1952 EEold(2,2,-i)=-b(10,i)+b(11,i)
1953 EEold(2,1,-i)=-b(12,i)+b(13,i)
1954 EEold(1,2,-i)=-b(12,i)-b(13,i)
1955 write(iout,*) "TU DOCHODZE"
1961 !c ee(2,1,i)=ee(1,2,i)
1966 "Coefficients of the cumulants (independent of valence angles)"
1967 do i=-nloctyp+1,nloctyp-1
1968 write (iout,*) 'Type ',onelet(iloctyp(i))
1970 write(iout,'(2f10.5)') B(3,i),B(5,i)
1972 write(iout,'(2f10.5)') B(2,i),B(4,i)
1975 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1979 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1983 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1992 ! Read torsional parameters in old format
1994 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1996 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1997 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1998 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2000 !el from energy module--------
2001 allocate(v1(nterm_old,ntortyp,ntortyp))
2002 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2003 !el---------------------------
2008 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2014 write (iout,'(/a/)') 'Torsional constants:'
2017 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2018 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2024 ! Read torsional parameters
2026 IF (TOR_MODE.eq.0) THEN
2027 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2028 read (itorp,*,end=113,err=113) ntortyp
2029 !el from energy module---------
2030 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2031 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2033 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2034 allocate(vlor2(maxlor,ntortyp,ntortyp))
2035 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2036 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2038 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2039 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2040 !el---------------------------
2043 !el---------------------------
2045 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2047 itortyp(i)=-itortyp(-i)
2049 itortyp(ntyp1)=ntortyp
2050 itortyp(-ntyp1)=-ntortyp
2052 write (iout,*) 'ntortyp',ntortyp
2054 do j=-ntortyp+1,ntortyp-1
2055 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2057 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2058 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2061 do k=1,nterm(i,j,iblock)
2062 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2064 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2065 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2066 v0ij=v0ij+si*v1(k,i,j,iblock)
2068 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2069 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2070 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2072 do k=1,nlor(i,j,iblock)
2073 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2074 vlor2(k,i,j),vlor3(k,i,j)
2075 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2078 v0(-i,-j,iblock)=v0ij
2084 write (iout,'(/a/)') 'Torsional constants:'
2086 do i=-ntortyp,ntortyp
2087 do j=-ntortyp,ntortyp
2088 write (iout,*) 'ityp',i,' jtyp',j
2089 write (iout,*) 'Fourier constants'
2090 do k=1,nterm(i,j,iblock)
2091 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2094 write (iout,*) 'Lorenz constants'
2095 do k=1,nlor(i,j,iblock)
2096 write (iout,'(3(1pe15.5))') &
2097 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2103 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2105 ! 6/23/01 Read parameters for double torsionals
2107 !el from energy module------------
2108 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2109 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2110 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2112 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2113 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2114 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2115 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2116 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2117 !---------------------------------
2121 do j=-ntortyp+1,ntortyp-1
2122 do k=-ntortyp+1,ntortyp-1
2123 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2124 ! write (iout,*) "OK onelett",
2127 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2128 .or. t3.ne.toronelet(k)) then
2129 write (iout,*) "Error in double torsional parameter file",&
2132 call MPI_Finalize(Ierror)
2134 stop "Error in double torsional parameter file"
2136 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2137 ntermd_2(i,j,k,iblock)
2138 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2139 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2140 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2141 ntermd_1(i,j,k,iblock))
2142 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2143 ntermd_1(i,j,k,iblock))
2144 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2145 ntermd_1(i,j,k,iblock))
2146 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2147 ntermd_1(i,j,k,iblock))
2148 ! Martix of D parameters for one dimesional foureir series
2149 do l=1,ntermd_1(i,j,k,iblock)
2150 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2151 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2152 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2153 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2154 ! write(iout,*) "whcodze" ,
2155 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2157 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2158 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2159 v2s(m,l,i,j,k,iblock),&
2160 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2161 ! Martix of D parameters for two dimesional fourier series
2162 do l=1,ntermd_2(i,j,k,iblock)
2164 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2165 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2166 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2167 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2176 write (iout,*) 'Constants for double torsionals'
2179 do j=-ntortyp+1,ntortyp-1
2180 do k=-ntortyp+1,ntortyp-1
2181 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2182 ' nsingle',ntermd_1(i,j,k,iblock),&
2183 ' ndouble',ntermd_2(i,j,k,iblock)
2185 write (iout,*) 'Single angles:'
2186 do l=1,ntermd_1(i,j,k,iblock)
2187 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2188 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2189 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2190 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2193 write (iout,*) 'Pairs of angles:'
2194 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2195 do l=1,ntermd_2(i,j,k,iblock)
2196 write (iout,'(i5,20f10.5)') &
2197 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2200 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2201 do l=1,ntermd_2(i,j,k,iblock)
2202 write (iout,'(i5,20f10.5)') &
2203 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2204 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2214 itype2loc(i)=itortyp(i)
2218 ELSE IF (TOR_MODE.eq.1) THEN
2220 !C read valence-torsional parameters
2221 read (itorp,*,end=121,err=121) ntortyp
2223 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2225 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2226 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2229 itype2loc(i)=itortyp(i)
2233 itortyp(i)=-itortyp(-i)
2235 do i=-ntortyp+1,ntortyp-1
2236 do j=-ntortyp+1,ntortyp-1
2237 !C first we read the cos and sin gamma parameters
2238 read (itorp,'(13x,a)',end=121,err=121) string
2239 write (iout,*) i,j,string
2240 read (itorp,*,end=121,err=121) &
2241 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2242 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2243 do k=1,nterm_kcc(j,i)
2244 do l=1,nterm_kcc_Tb(j,i)
2245 do ll=1,nterm_kcc_Tb(j,i)
2246 read (itorp,*,end=121,err=121) ii,jj,kk, &
2247 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2255 !c AL 4/8/16: Calculate coefficients from one-body parameters
2257 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2258 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2259 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2260 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2261 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2264 print *,i,itortyp(i)
2265 itortyp(i)=itype2loc(i)
2268 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2270 do i=-ntortyp+1,ntortyp-1
2271 do j=-ntortyp+1,ntortyp-1
2274 do k=1,nterm_kcc_Tb(j,i)
2275 do l=1,nterm_kcc_Tb(j,i)
2276 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2277 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2278 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2279 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2282 do k=1,nterm_kcc_Tb(j,i)
2283 do l=1,nterm_kcc_Tb(j,i)
2285 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2286 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2287 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2288 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2290 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2291 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2292 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2293 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2297 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2301 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2306 if (tor_mode.gt.0 .and. lprint) then
2307 !c Print valence-torsional parameters
2308 write (iout,'(a)') &
2309 "Parameters of the valence-torsional potentials"
2310 do i=-ntortyp+1,ntortyp-1
2311 do j=-ntortyp+1,ntortyp-1
2312 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2313 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2314 do k=1,nterm_kcc(j,i)
2315 do l=1,nterm_kcc_Tb(j,i)
2316 do ll=1,nterm_kcc_Tb(j,i)
2317 write (iout,'(3i5,2f15.4)')&
2318 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2326 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2327 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2328 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2329 !el from energy module---------
2330 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2331 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2333 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2334 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2335 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2336 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2338 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2339 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2340 !el---------------------------
2343 !el--------------------
2344 read (itorp_nucl,*,end=113,err=113) &
2345 (itortyp_nucl(i),i=1,ntyp_molec(2))
2346 ! print *,itortyp_nucl(:)
2347 !c write (iout,*) 'ntortyp',ntortyp
2350 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2351 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2354 do k=1,nterm_nucl(i,j)
2355 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2356 v0ij=v0ij+si*v1_nucl(k,i,j)
2359 do k=1,nlor_nucl(i,j)
2360 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2361 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2362 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2368 ! Read of Side-chain backbone correlation parameters
2369 ! Modified 11 May 2012 by Adasko
2372 read (isccor,*,end=119,err=119) nsccortyp
2374 !el from module energy-------------
2375 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2376 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2377 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2378 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2379 !-----------------------------------
2381 !el from module energy-------------
2382 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2384 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2386 isccortyp(i)=-isccortyp(-i)
2388 iscprol=isccortyp(20)
2389 ! write (iout,*) 'ntortyp',ntortyp
2391 !c maxinter is maximum interaction sites
2392 !el from module energy---------
2393 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2394 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2395 -nsccortyp:nsccortyp))
2396 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2397 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2398 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2399 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2400 !-----------------------------------
2404 read (isccor,*,end=119,err=119) &
2405 nterm_sccor(i,j),nlor_sccor(i,j)
2411 nterm_sccor(-i,j)=nterm_sccor(i,j)
2412 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2413 nterm_sccor(i,-j)=nterm_sccor(i,j)
2414 do k=1,nterm_sccor(i,j)
2415 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2417 if (j.eq.iscprol) then
2418 if (i.eq.isccortyp(10)) then
2419 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2420 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2423 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2424 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2425 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2426 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2427 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2428 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2429 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2432 if (i.eq.isccortyp(10)) then
2433 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2434 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2436 if (j.eq.isccortyp(10)) then
2437 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2438 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2440 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2441 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2442 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2443 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2444 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2445 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2449 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2450 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2451 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2452 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2455 do k=1,nlor_sccor(i,j)
2456 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2457 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2458 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2459 (1+vlor3sccor(k,i,j)**2)
2461 v0sccor(l,i,j)=v0ijsccor
2462 v0sccor(l,-i,j)=v0ijsccor1
2463 v0sccor(l,i,-j)=v0ijsccor2
2464 v0sccor(l,-i,-j)=v0ijsccor3
2470 !el from module energy-------------
2471 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2473 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2474 ! write (iout,*) 'ntortyp',ntortyp
2476 !c maxinter is maximum interaction sites
2477 !el from module energy---------
2478 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2479 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2480 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2481 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2482 !-----------------------------------
2486 read (isccor,*,end=119,err=119) &
2487 nterm_sccor(i,j),nlor_sccor(i,j)
2491 do k=1,nterm_sccor(i,j)
2492 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2494 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2497 do k=1,nlor_sccor(i,j)
2498 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2499 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2500 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2501 (1+vlor3sccor(k,i,j)**2)
2503 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2512 write (iout,'(/a/)') 'Torsional constants:'
2515 write (iout,*) 'ityp',i,' jtyp',j
2516 write (iout,*) 'Fourier constants'
2517 do k=1,nterm_sccor(i,j)
2518 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2520 write (iout,*) 'Lorenz constants'
2521 do k=1,nlor_sccor(i,j)
2522 write (iout,'(3(1pe15.5))') &
2523 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2530 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2531 ! interaction energy of the Gly, Ala, and Pro prototypes.
2534 ! Read electrostatic-interaction parameters
2539 write (iout,'(/a)') 'Electrostatic interaction constants:'
2540 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2541 'IT','JT','APP','BPP','AEL6','AEL3'
2543 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2544 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2545 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2546 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2551 app (i,j)=epp(i,j)*rri*rri
2552 bpp (i,j)=-2.0D0*epp(i,j)*rri
2553 ael6(i,j)=elpp6(i,j)*4.2D0**6
2554 ael3(i,j)=elpp3(i,j)*4.2D0**3
2556 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2562 ! Read side-chain interaction parameters.
2564 !el from module energy - COMMON.INTERACT-------
2565 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2566 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2567 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2569 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2570 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2571 allocate(epslip(ntyp,ntyp))
2579 !--------------------------------
2581 read (isidep,*,end=117,err=117) ipot,expon
2582 if (ipot.lt.1 .or. ipot.gt.5) then
2583 write (iout,'(2a)') 'Error while reading SC interaction',&
2584 'potential file - unknown potential type.'
2586 call MPI_Finalize(Ierror)
2592 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2593 ', exponents are ',expon,2*expon
2594 ! goto (10,20,30,30,40) ipot
2596 !----------------------- LJ potential ---------------------------------
2598 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2599 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2600 (sigma0(i),i=1,ntyp)
2602 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2603 write (iout,'(a/)') 'The epsilon array:'
2604 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2605 write (iout,'(/a)') 'One-body parameters:'
2606 write (iout,'(a,4x,a)') 'residue','sigma'
2607 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2610 !----------------------- LJK potential --------------------------------
2612 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2613 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2614 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2616 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2617 write (iout,'(a/)') 'The epsilon array:'
2618 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2619 write (iout,'(/a)') 'One-body parameters:'
2620 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2621 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2625 !---------------------- GB or BP potential -----------------------------
2628 ! print *,"I AM in SCELE",scelemode
2629 if (scelemode.eq.0) then
2631 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2633 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2634 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2635 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2636 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2638 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2641 ! For the GB potential convert sigma'**2 into chi'
2644 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2648 write (iout,'(/a/)') 'Parameters of the BP potential:'
2649 write (iout,'(a/)') 'The epsilon array:'
2650 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2651 write (iout,'(/a)') 'One-body parameters:'
2652 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2654 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2655 chip(i),alp(i),i=1,ntyp)
2658 ! print *,ntyp,"NTYP"
2659 allocate(icharge(ntyp1))
2660 ! print *,ntyp,icharge(i)
2662 read (isidep,*) (icharge(i),i=1,ntyp)
2663 print *,ntyp,icharge(i)
2664 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2665 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2666 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2667 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2668 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2669 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2670 allocate(epsintab(ntyp,ntyp))
2671 allocate(dtail(2,ntyp,ntyp))
2672 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2673 allocate(wqdip(2,ntyp,ntyp))
2674 allocate(wstate(4,ntyp,ntyp))
2675 allocate(dhead(2,2,ntyp,ntyp))
2676 allocate(nstate(ntyp,ntyp))
2677 allocate(debaykap(ntyp,ntyp))
2679 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2680 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2684 ! write (*,*) "Im in ALAB", i, " ", j
2686 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2687 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2688 chis(i,j),chis(j,i), & !2 w tej linii
2689 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2690 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2691 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2692 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2693 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2694 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2695 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2696 IF ((LaTeX).and.(i.gt.24)) then
2697 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2698 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2699 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2700 chis(i,j),chis(j,i) !2 w tej linii
2702 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2706 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2710 IF ((LaTeX).and.(i.gt.24)) then
2711 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2712 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2713 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2714 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2715 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2716 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2717 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2724 sigma(i,j) = sigma(j,i)
2725 sigmap1(i,j)=sigmap1(j,i)
2726 sigmap2(i,j)=sigmap2(j,i)
2727 sigiso1(i,j)=sigiso1(j,i)
2728 sigiso2(i,j)=sigiso2(j,i)
2729 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2730 nstate(i,j) = nstate(j,i)
2731 dtail(1,i,j) = dtail(1,j,i)
2732 dtail(2,i,j) = dtail(2,j,i)
2734 alphasur(k,i,j) = alphasur(k,j,i)
2735 wstate(k,i,j) = wstate(k,j,i)
2736 alphiso(k,i,j) = alphiso(k,j,i)
2739 dhead(2,1,i,j) = dhead(1,1,j,i)
2740 dhead(2,2,i,j) = dhead(1,2,j,i)
2741 dhead(1,1,i,j) = dhead(2,1,j,i)
2742 dhead(1,2,i,j) = dhead(2,2,j,i)
2744 epshead(i,j) = epshead(j,i)
2745 sig0head(i,j) = sig0head(j,i)
2748 wqdip(k,i,j) = wqdip(k,j,i)
2751 wquad(i,j) = wquad(j,i)
2752 epsintab(i,j) = epsintab(j,i)
2753 debaykap(i,j)=debaykap(j,i)
2754 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2761 !--------------------- GBV potential -----------------------------------
2763 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2764 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2765 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2766 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2768 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2769 write (iout,'(a/)') 'The epsilon array:'
2770 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2771 write (iout,'(/a)') 'One-body parameters:'
2772 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2773 's||/s_|_^2',' chip ',' alph '
2774 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2775 sigii(i),chip(i),alp(i),i=1,ntyp)
2778 write(iout,*)"Wrong ipot"
2784 !-----------------------------------------------------------------------
2785 ! Calculate the "working" parameters of SC interactions.
2787 !el from module energy - COMMON.INTERACT-------
2788 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2789 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2790 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2791 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2792 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2793 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2795 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2801 if (scelemode.eq.0) then
2810 sc_aa_tube_par(:)=0.0d0
2811 sc_bb_tube_par(:)=0.0d0
2813 !--------------------------------
2818 epslip(i,j)=epslip(j,i)
2821 if (scelemode.eq.0) then
2824 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2825 sigma(j,i)=sigma(i,j)
2826 rs0(i,j)=dwa16*sigma(i,j)
2831 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2832 'Working parameters of the SC interactions:',&
2833 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2838 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2840 ! print *,"SIGMA ZLA?",sigma(i,j)
2848 sigeps=dsign(1.0D0,epsij)
2850 aa_aq(i,j)=epsij*rrij*rrij
2851 ! print *,"ADASKO",epsij,rrij,expon
2852 bb_aq(i,j)=-sigeps*epsij*rrij
2853 aa_aq(j,i)=aa_aq(i,j)
2854 bb_aq(j,i)=bb_aq(i,j)
2855 epsijlip=epslip(i,j)
2856 sigeps=dsign(1.0D0,epsijlip)
2857 epsijlip=dabs(epsijlip)
2858 aa_lip(i,j)=epsijlip*rrij*rrij
2859 bb_lip(i,j)=-sigeps*epsijlip*rrij
2860 aa_lip(j,i)=aa_lip(i,j)
2861 bb_lip(j,i)=bb_lip(i,j)
2862 !C write(iout,*) aa_lip
2863 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2864 sigt1sq=sigma0(i)**2
2865 sigt2sq=sigma0(j)**2
2868 ratsig1=sigt2sq/sigt1sq
2869 ratsig2=1.0D0/ratsig1
2870 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2871 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2872 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2876 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2877 sigmaii(i,j)=rsum_max
2878 sigmaii(j,i)=rsum_max
2880 ! sigmaii(i,j)=r0(i,j)
2881 ! sigmaii(j,i)=r0(i,j)
2883 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2884 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2885 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2886 augm(i,j)=epsij*r_augm**(2*expon)
2887 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2894 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2895 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2896 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2901 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2902 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2903 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2904 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2905 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2906 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2907 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2908 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2909 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2910 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2911 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2912 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2913 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2914 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2915 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2916 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2925 read (isidep_nucl,*) ipot_nucl
2926 ! print *,"TU?!",ipot_nucl
2927 if (ipot_nucl.eq.1) then
2928 do i=1,ntyp_molec(2)
2929 do j=i,ntyp_molec(2)
2930 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2931 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2935 do i=1,ntyp_molec(2)
2936 do j=i,ntyp_molec(2)
2937 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2938 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2939 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2943 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2944 do i=1,ntyp_molec(2)
2945 do j=i,ntyp_molec(2)
2946 rrij=sigma_nucl(i,j)
2950 epsij=4*eps_nucl(i,j)
2951 sigeps=dsign(1.0D0,epsij)
2953 aa_nucl(i,j)=epsij*rrij*rrij
2954 bb_nucl(i,j)=-sigeps*epsij*rrij
2955 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2956 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2957 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2958 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2959 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2960 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2963 aa_nucl(i,j)=aa_nucl(j,i)
2964 bb_nucl(i,j)=bb_nucl(j,i)
2965 ael3_nucl(i,j)=ael3_nucl(j,i)
2966 ael6_nucl(i,j)=ael6_nucl(j,i)
2967 ael63_nucl(i,j)=ael63_nucl(j,i)
2968 ael32_nucl(i,j)=ael32_nucl(j,i)
2969 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2970 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2971 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2972 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2973 eps_nucl(i,j)=eps_nucl(j,i)
2974 sigma_nucl(i,j)=sigma_nucl(j,i)
2975 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2979 write(iout,*) "tube param"
2980 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2981 ccavtubpep,dcavtubpep,tubetranenepep
2982 sigmapeptube=sigmapeptube**6
2983 sigeps=dsign(1.0D0,epspeptube)
2984 epspeptube=dabs(epspeptube)
2985 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2986 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2987 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2989 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2990 ccavtub(i),dcavtub(i),tubetranene(i)
2991 sigmasctube=sigmasctube**6
2992 sigeps=dsign(1.0D0,epssctube)
2993 epssctube=dabs(epssctube)
2994 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2995 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2996 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2998 !-----------------READING SC BASE POTENTIALS-----------------------------
2999 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3000 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3001 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3002 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3003 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3004 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3005 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3006 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3007 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3008 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3009 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3010 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3011 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3012 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3013 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3014 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3016 write (iout,*) "ESCBASEPARM"
3017 do i=1,ntyp_molec(1)
3018 do j=1,ntyp_molec(2) ! without U then we will take T for U
3019 ! write (*,*) "Im in ", i, " ", j
3020 read(isidep_scbase,*) &
3021 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3022 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3023 ! write(*,*) "eps",eps_scbase(i,j)
3024 read(isidep_scbase,*) &
3025 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3026 chis_scbase(i,j,1),chis_scbase(i,j,2)
3027 read(isidep_scbase,*) &
3028 dhead_scbasei(i,j), &
3029 dhead_scbasej(i,j), &
3030 rborn_scbasei(i,j),rborn_scbasej(i,j)
3031 read(isidep_scbase,*) &
3032 (wdipdip_scbase(k,i,j),k=1,3), &
3033 (wqdip_scbase(k,i,j),k=1,2)
3034 read(isidep_scbase,*) &
3035 alphapol_scbase(i,j), &
3036 epsintab_scbase(i,j)
3037 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3038 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3039 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3040 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3041 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3042 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3043 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3044 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3046 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3047 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3048 write(*,*) "eps",eps_scbase(i,j)
3050 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3051 chis_scbase(i,j,1),chis_scbase(i,j,2)
3053 dhead_scbasei(i,j), &
3054 dhead_scbasej(i,j), &
3055 rborn_scbasei(i,j),rborn_scbasej(i,j)
3057 (wdipdip_scbase(k,i,j),k=1,3), &
3058 (wqdip_scbase(k,i,j),k=1,2)
3060 alphapol_scbase(i,j), &
3061 epsintab_scbase(i,j)
3066 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3067 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3068 write(*,*) "eps",eps_scbase(i,j)
3070 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3071 chis_scbase(i,j,1),chis_scbase(i,j,2)
3073 dhead_scbasei(i,j), &
3074 dhead_scbasej(i,j), &
3075 rborn_scbasei(i,j),rborn_scbasej(i,j)
3077 (wdipdip_scbase(k,i,j),k=1,3), &
3078 (wqdip_scbase(k,i,j),k=1,2)
3080 alphapol_scbase(i,j), &
3081 epsintab_scbase(i,j)
3084 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3085 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3087 do i=1,ntyp_molec(1)
3088 do j=1,ntyp_molec(2)-1
3089 epsij=eps_scbase(i,j)
3090 rrij=sigma_scbase(i,j)
3095 sigeps=dsign(1.0D0,epsij)
3097 aa_scbase(i,j)=epsij*rrij*rrij
3098 bb_scbase(i,j)=-sigeps*epsij*rrij
3101 !-----------------READING PEP BASE POTENTIALS-------------------
3102 allocate(eps_pepbase(ntyp_molec(2)))
3103 allocate(sigma_pepbase(ntyp_molec(2)))
3104 allocate(chi_pepbase(ntyp_molec(2),2))
3105 allocate(chipp_pepbase(ntyp_molec(2),2))
3106 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3107 allocate(sigmap1_pepbase(ntyp_molec(2)))
3108 allocate(sigmap2_pepbase(ntyp_molec(2)))
3109 allocate(chis_pepbase(ntyp_molec(2),2))
3110 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3112 write (iout,*) "EPEPBASEPARM"
3114 do j=1,ntyp_molec(2) ! without U then we will take T for U
3115 write (*,*) "Im in ", i, " ", j
3116 read(isidep_pepbase,*) &
3117 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3118 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3119 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3120 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3121 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3122 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3123 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3124 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3125 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3126 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3128 write(*,*) "eps",eps_pepbase(j)
3129 read(isidep_pepbase,*) &
3130 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3131 chis_pepbase(j,1),chis_pepbase(j,2)
3132 read(isidep_pepbase,*) &
3133 (wdipdip_pepbase(k,j),k=1,3)
3134 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3135 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3137 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3138 chis_pepbase(j,1),chis_pepbase(j,2)
3140 (wdipdip_pepbase(k,j),k=1,3)
3144 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3145 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3147 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3148 chis_pepbase(j,1),chis_pepbase(j,2)
3150 (wdipdip_pepbase(k,j),k=1,3)
3152 allocate(aa_pepbase(ntyp_molec(2)))
3153 allocate(bb_pepbase(ntyp_molec(2)))
3155 do j=1,ntyp_molec(2)-1
3156 epsij=eps_pepbase(j)
3157 rrij=sigma_pepbase(j)
3162 sigeps=dsign(1.0D0,epsij)
3164 aa_pepbase(j)=epsij*rrij*rrij
3165 bb_pepbase(j)=-sigeps*epsij*rrij
3167 !--------------READING SC PHOSPHATE-------------------------------------
3168 allocate(eps_scpho(ntyp_molec(1)))
3169 allocate(sigma_scpho(ntyp_molec(1)))
3170 allocate(chi_scpho(ntyp_molec(1),2))
3171 allocate(chipp_scpho(ntyp_molec(1),2))
3172 allocate(alphasur_scpho(4,ntyp_molec(1)))
3173 allocate(sigmap1_scpho(ntyp_molec(1)))
3174 allocate(sigmap2_scpho(ntyp_molec(1)))
3175 allocate(chis_scpho(ntyp_molec(1),2))
3176 allocate(wqq_scpho(ntyp_molec(1)))
3177 allocate(wqdip_scpho(2,ntyp_molec(1)))
3178 allocate(alphapol_scpho(ntyp_molec(1)))
3179 allocate(epsintab_scpho(ntyp_molec(1)))
3180 allocate(dhead_scphoi(ntyp_molec(1)))
3181 allocate(rborn_scphoi(ntyp_molec(1)))
3182 allocate(rborn_scphoj(ntyp_molec(1)))
3183 allocate(alphi_scpho(ntyp_molec(1)))
3187 do j=1,ntyp_molec(1) ! without U then we will take T for U
3188 write (*,*) "Im in scpho ", i, " ", j
3189 read(isidep_scpho,*) &
3190 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3191 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3192 write(*,*) "eps",eps_scpho(j)
3193 read(isidep_scpho,*) &
3194 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3195 chis_scpho(j,1),chis_scpho(j,2)
3196 read(isidep_scpho,*) &
3197 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3198 read(isidep_scpho,*) &
3199 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3201 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3202 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3203 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3204 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3205 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3206 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3207 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3208 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3212 allocate(aa_scpho(ntyp_molec(1)))
3213 allocate(bb_scpho(ntyp_molec(1)))
3215 do j=1,ntyp_molec(1)
3222 sigeps=dsign(1.0D0,epsij)
3224 aa_scpho(j)=epsij*rrij*rrij
3225 bb_scpho(j)=-sigeps*epsij*rrij
3229 read(isidep_peppho,*) &
3230 eps_peppho,sigma_peppho
3231 read(isidep_peppho,*) &
3232 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3233 read(isidep_peppho,*) &
3234 (wqdip_peppho(k),k=1,2)
3242 sigeps=dsign(1.0D0,epsij)
3244 aa_peppho=epsij*rrij*rrij
3245 bb_peppho=-sigeps*epsij*rrij
3248 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3253 ! Define the SC-p interaction constants (hard-coded; old style)
3256 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3258 ! aad(i,1)=0.3D0*4.0D0**12
3259 ! Following line for constants currently implemented
3260 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3261 aad(i,1)=1.5D0*4.0D0**12
3262 ! aad(i,1)=0.17D0*5.6D0**12
3264 ! "Soft" SC-p repulsion
3266 ! Following line for constants currently implemented
3267 ! aad(i,1)=0.3D0*4.0D0**6
3268 ! "Hard" SC-p repulsion
3269 bad(i,1)=3.0D0*4.0D0**6
3270 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3279 ! 8/9/01 Read the SC-p interaction constants from file
3282 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3285 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3286 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3287 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3288 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3292 write (iout,*) "Parameters of SC-p interactions:"
3294 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3295 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3300 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3302 do i=1,ntyp_molec(2)
3303 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3305 do i=1,ntyp_molec(2)
3306 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3307 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3309 r0pp=1.12246204830937298142*5.16158
3315 ! Define the constants of the disulfide bridge
3319 ! Old arbitrary potential - commented out.
3324 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3325 ! energy surface of diethyl disulfide.
3326 ! A. Liwo and U. Kozlowska, 11/24/03
3342 allocate(ichargelipid(ntyp_molec(4)))
3343 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3344 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3345 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3346 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3347 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3348 kjtokcal=0.2390057361
3350 !HERE THE MASS of MARTINI
3351 write(*,*) "before MARTINI PARAM"
3352 do i=1,ntyp_molec(4)
3358 !relative dielectric constant = 15 for implicit screening
3359 k_coulomb_lip=332.0d0/15.0d0
3360 !kbond = 1250 kJ/(mol*nm*2)
3361 kbondlip=1250.0d0*kjtokcal/100.0d0
3363 lip_angle_force=0.0d0
3364 if (DRY_MARTINI.gt.0) then
3365 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3366 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3367 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3368 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3369 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3370 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3371 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3372 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3374 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3375 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3376 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3377 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3378 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3379 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3380 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3381 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3383 lip_angle_angle=0.0d0
3384 lip_angle_angle(3,12,12)=120.0/krad
3385 lip_angle_angle(3,12,18)=180.0/krad
3386 lip_angle_angle(3,18,16)=180.0/krad
3387 lip_angle_angle(12,18,16)=180.0/krad
3388 lip_angle_angle(18,16,18)=120.0/krad
3389 lip_angle_angle(16,18,18)=180.0/krad
3390 lip_angle_angle(12,18,18)=180.0/krad
3391 lip_angle_angle(18,18,18)=180.0/krad
3392 read(ilipbond,*) temp1
3394 read(ilipbond,*) temp1, lip_bond(i,1), &
3395 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3396 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3397 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3398 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3401 lip_bond(i,j)=lip_bond(i,j)*10
3405 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3406 read(ilipnonbond,*) temp1
3408 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3409 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3410 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3411 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3412 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3414 ! write(*,*) i, lip_eps(i,18)
3416 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3419 read(ilipnonbond,*) temp1
3421 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3422 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3423 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3424 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3425 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3428 lip_sig(i,j)=lip_sig(i,j)*10.0
3431 write(*,*) "after MARTINI PARAM"
3435 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3436 allocate(alphapolcat2(ntyp,ntyp))
3437 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3438 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3439 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3440 allocate(epsintabcat(ntyp,ntyp))
3441 allocate(dtailcat(2,ntyp,ntyp))
3442 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3443 allocate(wqdipcat(2,ntyp,ntyp))
3444 allocate(wstatecat(4,ntyp,ntyp))
3445 allocate(dheadcat(2,2,ntyp,ntyp))
3446 allocate(nstatecat(ntyp,ntyp))
3447 allocate(debaykapcat(ntyp,ntyp))
3449 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3450 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3451 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3452 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3453 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3456 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3457 write(*,*) "before ions",oldion
3459 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3460 if (oldion.eq.0) then
3461 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3462 allocate(icharge(1:ntyp1))
3463 read(iion,*) (icharge(i),i=1,ntyp)
3467 print *,ntyp_molec(5)
3468 do i=1,ntyp_molec(5)
3469 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3470 print *,msc(i,5),restok(i,5)
3476 do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3478 ! do j=1,ntyp_molec(5)
3479 ! write (*,*) "Im in ALAB", i, " ", j
3481 epscat(i,j),sigmacat(i,j), &
3482 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3483 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3485 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3486 ! chiscat(i,j),chiscat(j,i), &
3487 chis1cat(i,j),chis2cat(i,j), &
3489 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3490 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3491 dtailcat(1,i,j),dtailcat(2,i,j), &
3492 epsheadcat(i,j),sig0headcat(i,j), &
3494 ! rborncat(i,j),rborncat(j,i),&
3495 rborn1cat(i,j),rborn2cat(i,j),&
3496 (wqdipcat(k,i,j),k=1,2), &
3497 alphapolcat(i,j),alphapolcat2(j,i), &
3498 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3500 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3501 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3503 ! write (iout,*) 'i= ', i, ' j= ', j
3504 ! write (iout,*) 'epsi0= ', epscat(1,j)
3505 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3506 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3507 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3508 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3509 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3510 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3511 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3512 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3513 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3514 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3515 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3516 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3517 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3518 ! write (iout,*) 'a1= ', rborncat(j,1)
3519 ! write (iout,*) 'a2= ', rborncat(1,j)
3520 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3521 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3522 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3523 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3527 ! If ((i.eq.1).and.(j.eq.27)) then
3528 ! write (iout,*) 'SEP'
3529 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3530 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3535 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3537 do j=1,ntyp_molec(5)
3541 sigeps=dsign(1.0D0,epsij)
3543 aa_aq_cat(i,j)=epsij*rrij*rrij
3544 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3549 do j=1,ntyp_molec(5)-1
3551 write (iout,*) 'i= ', i, ' j= ', j
3552 write (iout,*) 'epsi0= ', epscat(i,j)
3553 write (iout,*) 'sigma0= ', sigmacat(i,j)
3554 write (iout,*) 'chi1= ', chi1cat(i,j)
3555 write (iout,*) 'chi1= ', chi2cat(i,j)
3556 write (iout,*) 'chip1= ', chipp1cat(i,j)
3557 write (iout,*) 'chip2= ', chipp2cat(i,j)
3558 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3559 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3560 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3561 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3562 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3563 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3564 write (iout,*) 'chis1= ', chis1cat(i,j)
3565 write (iout,*) 'chis1= ', chis2cat(i,j)
3566 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3567 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3568 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3569 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3570 write (iout,*) 'a1= ', rborn1cat(i,j)
3571 write (iout,*) 'a2= ', rborn2cat(i,j)
3572 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3573 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3574 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3575 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3576 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3577 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3580 If ((i.eq.1).and.(j.eq.27)) then
3581 write (iout,*) 'SEP'
3582 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3583 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3590 ! read number of Zn2+
3591 ! here two denotes the Zn2+ and Cu2+
3592 write(iout,*) "before TRANPARM"
3593 allocate(aomicattr(0:3,2))
3594 allocate(athetacattran(0:6,5,2))
3595 allocate(agamacattran(3,5,2))
3596 allocate(acatshiftdsc(5,2))
3597 allocate(bcatshiftdsc(5,2))
3598 allocate(demorsecat(5,2))
3599 allocate(alphamorsecat(5,2))
3600 allocate(x0catleft(5,2))
3601 allocate(x0catright(5,2))
3602 allocate(x0cattrans(5,2))
3603 allocate(ntrantyp(2))
3604 do i=1,1 ! currently only Zn2+
3606 read(iiontran,*) ntrantyp(i)
3609 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3610 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3611 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3612 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3613 read(iiontran,*) (aomicattr(j,i),j=0,3)
3615 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3616 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3617 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3622 write (iout,'(/a)') "Disulfide bridge parameters:"
3623 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3624 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3625 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3626 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3629 if (shield_mode.gt.0) then
3630 pi=4.0D0*datan(1.0D0)
3631 !C VSolvSphere the volume of solving sphere
3633 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3634 !C there will be no distinction between proline peptide group and normal peptide
3635 !C group in case of shielding parameters
3636 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3637 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3638 write (iout,*) VSolvSphere,VSolvSphere_div
3639 !C long axis of side chain
3641 long_r_sidechain(i)=vbldsc0(1,i)
3642 ! if (scelemode.eq.0) then
3643 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3644 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3646 ! short_r_sidechain(i)=sigma(i,i)
3648 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3655 111 write (iout,*) "Error reading bending energy parameters."
3657 112 write (iout,*) "Error reading rotamer energy parameters."
3659 113 write (iout,*) "Error reading torsional energy parameters."
3661 114 write (iout,*) "Error reading double torsional energy parameters."
3663 115 write (iout,*) &
3664 "Error reading cumulant (multibody energy) parameters."
3666 116 write (iout,*) "Error reading electrostatic energy parameters."
3668 117 write (iout,*) "Error reading side chain interaction parameters."
3670 118 write (iout,*) "Error reading SCp interaction parameters."
3672 119 write (iout,*) "Error reading SCCOR parameters"
3674 121 write (iout,*) "Error in Czybyshev parameters"
3676 123 write(iout,*) "Error in transition metal parameters"
3679 call MPI_Finalize(Ierror)
3683 end subroutine parmread
3685 !-----------------------------------------------------------------------------
3687 !-----------------------------------------------------------------------------
3688 subroutine printmat(ldim,m,n,iout,key,a)
3691 character(len=3),dimension(n) :: key
3692 real(kind=8),dimension(ldim,n) :: a
3694 integer :: i,j,k,m,iout,nlim
3698 write (iout,1000) (key(k),k=i,nlim)
3700 1000 format (/5x,8(6x,a3))
3701 1020 format (/80(1h-)/)
3703 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3706 1010 format (a3,2x,8(f9.4))
3708 end subroutine printmat
3709 !-----------------------------------------------------------------------------
3711 !-----------------------------------------------------------------------------
3713 ! Read the PDB file and convert the peptide geometry into virtual-chain
3716 use energy_data, only: itype,istype
3720 ! use control, only: rescode,sugarcode
3721 ! implicit real*8 (a-h,o-z)
3722 ! include 'DIMENSIONS'
3723 ! include 'COMMON.LOCAL'
3724 ! include 'COMMON.VAR'
3725 ! include 'COMMON.CHAIN'
3726 ! include 'COMMON.INTERACT'
3727 ! include 'COMMON.IOUNITS'
3728 ! include 'COMMON.GEO'
3729 ! include 'COMMON.NAMES'
3730 ! include 'COMMON.CONTROL'
3731 ! include 'COMMON.DISTFIT'
3732 ! include 'COMMON.SETUP'
3733 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3735 logical :: lprn=.true.,fail
3736 real(kind=8),dimension(3) :: e1,e2,e3
3737 real(kind=8) :: dcj,efree_temp
3738 character(len=3) :: seq,res,res2
3739 character(len=5) :: atom
3740 character(len=80) :: card
3741 real(kind=8),dimension(3,20) :: sccor
3742 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3743 integer :: isugar,molecprev,firstion
3744 character*1 :: sugar
3746 real(kind=8),dimension(3) :: ccc
3748 integer,dimension(2,maxres/3) :: hfrag_alloc
3749 integer,dimension(4,maxres/3) :: bfrag_alloc
3750 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3751 real(kind=8),dimension(:,:), allocatable :: c_temporary
3752 integer,dimension(:,:) , allocatable :: itype_temporary
3753 integer,dimension(:),allocatable :: istype_temp
3760 ! write (2,*) "UNRES_PDB",unres_pdb
3780 !-----------------------------
3781 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3782 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3783 if(.not. allocated(istype)) allocate(istype(maxres))
3785 read (ipdbin,'(a80)',end=10) card
3786 write (iout,'(a)') card
3787 if (card(:5).eq.'HELIX') then
3790 read(card(22:25),*) hfrag(1,nhfrag)
3791 read(card(34:37),*) hfrag(2,nhfrag)
3793 if (card(:5).eq.'SHEET') then
3796 read(card(24:26),*) bfrag(1,nbfrag)
3797 read(card(35:37),*) bfrag(2,nbfrag)
3798 !rc----------------------------------------
3799 !rc to be corrected !!!
3800 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3801 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3802 !rc----------------------------------------
3804 if (card(:3).eq.'END') then
3806 else if (card(:3).eq.'TER') then
3811 itype(ires_old,molecule)=ntyp1_molec(molecule)
3812 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3813 nres_molec(molecule)=nres_molec(molecule)+2
3814 ! if (molecule.eq.4) ires=ires+2
3816 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3819 dc(j,ires)=sccor(j,iii)
3822 call sccenter(ires,iii,sccor)
3826 else if (card(:3).eq.'BRA') then
3830 itype(ires,molecule)=ntyp1_molec(molecule)-1
3831 nres_molec(molecule)=nres_molec(molecule)+1
3835 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3836 ! Fish out the ATOM cards.
3837 ! write(iout,*) 'card',card(1:20)
3838 ! print *,"ATU ",card(1:6), CARD(3:6)
3840 if (index(card(1:4),'ATOM').gt.0) then
3841 read (card(12:16),*) atom
3842 ! write (iout,*) "! ",atom," !",ires
3843 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3844 if (card(14:16).eq.'LIP') then
3849 nres_molec(molecule)=nres_molec(molecule)+1
3850 itype(ires,molecule)=ntyp1_molec(molecule)
3859 nres_molec(molecule)=nres_molec(molecule)+1
3860 read (card(18:20),'(a3)') res
3862 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3863 if (UNRES_PDB) then!
3864 itype(ires,molecule)=rescode(ires,res,0,molecule)
3866 itype(ires,molecule)=rescode_lip(res,ires)
3869 read (card(23:26),*) ires
3870 read (card(18:20),'(a3)') res
3871 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3872 ! & " ires_old",ires_old
3873 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3874 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3875 if (ires-ishift+ishift1.ne.ires_old) then
3876 ! Calculate the CM of the preceding residue.
3877 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3879 ! write (iout,*) "Calculating sidechain center iii",iii
3882 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3885 call sccenter(ires_old,iii,sccor)
3889 ! Start new residue.
3890 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3893 else if (ibeg.eq.1) then
3894 write (iout,*) "BEG ires",ires
3896 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3899 nres_molec(molecule)=nres_molec(molecule)+1
3901 ires=ires-ishift+ishift1
3903 ! write (iout,*) "ishift",ishift," ires",ires,&
3904 ! " ires_old",ires_old
3906 else if (ibeg.eq.2) then
3908 ishift=-ires_old+ires-1 !!!!!
3909 ishift1=ishift1-1 !!!!!
3910 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3911 ires=ires-ishift+ishift1
3912 ! print *,ires,ishift,ishift1
3916 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3917 ires=ires-ishift+ishift1
3920 ! print *,'atom',ires,atom
3921 if (res.eq.'ACE' .or. res.eq.'NHE') then
3924 if (atom.eq.'CA '.or.atom.eq.'N ') then
3926 itype(ires,molecule)=rescode(ires,res,0,molecule)
3928 ! nres_molec(molecule)=nres_molec(molecule)+1
3932 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3933 ! nres_molec(molecule)=nres_molec(molecule)+1
3934 read (card(19:19),'(a1)') sugar
3935 isugar=sugarcode(sugar,ires)
3936 ! if (ibeg.eq.1) then
3940 ! print *,"ires=",ires,istype(ires)
3946 ires=ires-ishift+ishift1
3948 ! write (iout,*) "ires_old",ires_old," ires",ires
3949 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3952 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3953 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3954 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3955 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3956 ! print *,ires,ishift,ishift1
3957 ! write (iout,*) "backbone ",atom
3959 write (iout,'(2i3,2x,a,3f8.3)') &
3960 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3963 nres_molec(molecule)=nres_molec(molecule)+1
3965 sccor(j,iii)=c(j,ires)
3967 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3968 atom.eq."C2'" .or. atom.eq."C3'" &
3969 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3970 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3971 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3972 ! print *,ires,ishift,ishift1
3976 ! sccor(j,iii)=c(j,ires)
3979 c(j,ires)=c(j,ires)+ccc(j)/5.0
3981 print *,counter,molecule
3982 if (counter.eq.5) then
3984 nres_molec(molecule)=nres_molec(molecule)+1
3987 ! sccor(j,iii)=c(j,ires)
3991 ! print *, "ATOM",atom(1:3)
3992 else if (atom.eq."C5'") then
3993 read (card(19:19),'(a1)') sugar
3994 isugar=sugarcode(sugar,ires)
3999 ! print *,ires,istype(ires)
4003 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4004 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4005 nres_molec(molecule)=nres_molec(molecule)+1
4006 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4010 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4012 else if ((atom.eq."C1'").and.unres_pdb) then
4014 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4015 ! write (*,*) card(23:27),ires,itype(ires,1)
4016 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4017 atom.ne.'N' .and. atom.ne.'C' .and. &
4018 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4019 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4020 .and. atom.ne.'P '.and. &
4021 atom(1:1).ne.'H' .and. &
4022 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4024 ! write (iout,*) "sidechain ",atom
4025 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4026 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4027 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4029 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4033 ! print *,"IONS",ions,card(1:6)
4034 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4035 if (firstion.eq.0) then
4039 dc(j,ires)=sccor(j,iii)
4042 call sccenter(ires,iii,sccor)
4045 read (card(12:16),*) atom
4046 ! print *,"HETATOM", atom
4047 read (card(18:20),'(a3)') res
4048 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4049 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4050 .or.(atom(1:2).eq.'K '.or.(atom(1:2).eq.'ZN'))) &
4053 if (molecule.ne.5) molecprev=molecule
4055 nres_molec(molecule)=nres_molec(molecule)+1
4056 print *,"HERE",nres_molec(molecule)
4058 itype(ires,molecule)=rescode(ires,res,0,molecule)
4059 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4063 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4064 if (ires.eq.0) return
4065 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4068 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4070 nres_molec(molecule)=nres_molec(molecule)-2
4071 print *,'I have',nres, nres_molec(:)
4073 do k=1,4 ! ions are without dummy
4074 if (nres_molec(k).eq.0) cycle
4076 ! write (iout,*) i,itype(i,1)
4077 ! if (itype(i,1).eq.ntyp1) then
4078 ! write (iout,*) "dummy",i,itype(i,1)
4080 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4081 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4085 if (itype(i,k).eq.ntyp1_molec(k)) then
4086 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4087 if (itype(i+2,k).eq.0) then
4088 ! print *,"masz sieczke"
4090 if (itype(i+2,j).ne.0) then
4092 itype(i+1,j)=ntyp1_molec(j)
4093 nres_molec(k)=nres_molec(k)-1
4094 nres_molec(j)=nres_molec(j)+1
4100 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4101 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4102 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4103 ! if (unres_pdb) then
4104 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4105 ! print *,i,'tu dochodze'
4106 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4114 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4118 dcj=(c(j,i-2)-c(j,i-3))/2.0
4119 if (dcj.eq.0) dcj=1.23591524223
4124 else !itype(i+1,1).eq.ntyp1
4125 ! if (unres_pdb) then
4126 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4127 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4134 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4135 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4139 dcj=(c(j,i+3)-c(j,i+2))/2.0
4140 if (dcj.eq.0) dcj=1.23591524223
4145 endif !itype(i+1,1).eq.ntyp1
4146 endif !itype.eq.ntyp1
4150 ! Calculate the CM of the last side chain.
4154 dc(j,ires)=sccor(j,iii)
4157 call sccenter(ires,iii,sccor)
4163 ! print *,"molecule",molecule
4164 if ((itype(nres,1).ne.10)) then
4167 if (molecule.eq.5) molecule=molecprev
4168 itype(nres,molecule)=ntyp1_molec(molecule)
4169 nres_molec(molecule)=nres_molec(molecule)+1
4170 ! if (unres_pdb) then
4171 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4172 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4179 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4183 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4184 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4185 c(j,2*nres)=c(j,nres)
4189 ! print *,'I have',nres, nres_molec(:)
4191 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4193 if (nres.ne.nres0) then
4194 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4196 stop "Error nres value in WHAM input"
4199 !---------------------------------
4200 !el reallocate tables
4203 ! hfrag_alloc(j,i)=hfrag(j,i)
4206 ! bfrag_alloc(j,i)=bfrag(j,i)
4212 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4213 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4214 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4215 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4219 ! hfrag(j,i)=hfrag_alloc(j,i)
4224 ! bfrag(j,i)=bfrag_alloc(j,i)
4227 !el end reallocate tables
4228 !---------------------------------
4236 c(j,2*nres)=c(j,nres)
4239 if (itype(1,1).eq.ntyp1) then
4243 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4244 call refsys(2,3,4,e1,e2,e3,fail)
4251 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4252 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4256 dcj=(c(j,4)-c(j,3))/2.0
4262 ! First lets assign correct dummy to correct type of chain
4264 if (itype(1,1).eq.ntyp1) then
4265 if (itype(2,1).eq.0) then
4267 if (itype(2,j).ne.0) then
4269 itype(1,j)=ntyp1_molec(j)
4270 nres_molec(1)=nres_molec(1)-1
4271 nres_molec(j)=nres_molec(j)+1
4278 print *,'I have',nres, nres_molec(:)
4280 ! Copy the coordinates to reference coordinates
4286 ! Calculate internal coordinates.
4288 write (iout,'(/a)') &
4289 "Cartesian coordinates of the reference structure"
4290 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4291 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4293 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4294 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4295 (c(j,ires+nres),j=1,3)
4298 ! znamy już nres wiec mozna alokowac tablice
4299 ! Calculate internal coordinates.
4300 if(me.eq.king.or..not.out1file)then
4301 write (iout,'(a)') &
4302 "Backbone and SC coordinates as read from the PDB"
4304 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4305 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4306 (c(j,nres+ires),j=1,3)
4309 ! NOW LETS ROCK! SORTING
4310 allocate(c_temporary(3,2*nres))
4311 allocate(itype_temporary(nres,5))
4312 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4313 if (.not.allocated(istype)) write(iout,*) &
4314 "SOMETHING WRONG WITH ISTYTPE"
4315 allocate(istype_temp(nres))
4316 itype_temporary(:,:)=0
4320 if (itype(i,k).ne.0) then
4322 c_temporary(j,seqalingbegin)=c(j,i)
4323 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4325 itype_temporary(seqalingbegin,k)=itype(i,k)
4326 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4327 istype_temp(seqalingbegin)=istype(i)
4328 molnum(seqalingbegin)=k
4329 seqalingbegin=seqalingbegin+1
4335 c(j,i)=c_temporary(j,i)
4337 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4343 itype(i,k)=itype_temporary(i,k)
4344 istype(i)=istype_temp(i)
4347 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4348 ! I have only ions now dummy atoms in the system
4350 itype(1,5)=itype(2,5)
4353 itype(i,5)=itype(i+1,5)
4357 nres_molec(1)=nres_molec(1)-1
4359 ! if (itype(1,1).eq.ntyp1) then
4362 ! if (unres_pdb) then
4363 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4364 ! call refsys(2,3,4,e1,e2,e3,fail)
4371 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4375 ! dcj=(c(j,4)-c(j,3))/2.0
4377 ! c(j,nres+1)=c(j,1)
4383 write (iout,'(/a)') &
4384 "Cartesian coordinates of the reference structure after sorting"
4385 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4386 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4388 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4389 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4390 (c(j,ires+nres),j=1,3)
4394 ! print *,seqalingbegin,nres
4395 if(.not.allocated(vbld)) then
4396 allocate(vbld(2*nres))
4401 if(.not.allocated(vbld_inv)) then
4402 allocate(vbld_inv(2*nres))
4408 if(.not.allocated(theta)) then
4409 allocate(theta(nres+2))
4413 if(.not.allocated(phi)) allocate(phi(nres+2))
4414 if(.not.allocated(alph)) allocate(alph(nres+2))
4415 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4416 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4417 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4418 if(.not.allocated(costtab)) allocate(costtab(nres))
4419 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4420 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4421 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4422 if(.not.allocated(xxref)) allocate(xxref(nres))
4423 if(.not.allocated(yyref)) allocate(yyref(nres))
4424 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4425 if(.not.allocated(dc_norm)) then
4426 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4427 allocate(dc_norm(3,0:2*nres+2))
4431 call int_from_cart(.true.,.false.)
4432 call sc_loc_geom(.false.)
4434 thetaref(i)=theta(i)
4444 dc(j,i)=c(j,i+1)-c(j,i)
4445 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4450 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4451 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4453 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4457 ! Copy the coordinates to reference coordinates
4458 ! Splits to single chain if occurs
4462 ! cref(j,i,cou)=c(j,i)
4466 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4467 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4468 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4469 !-----------------------------
4473 write (iout,*) "symetr", symetr
4476 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4478 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4481 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4486 cref(j,i,cou)=c(j,i)
4487 cref(j,i+nres,cou)=c(j,i+nres)
4489 chain_rep(j,lll,kkk)=c(j,i)
4490 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4494 write (iout,*) chain_length
4495 if (chain_length.eq.0) chain_length=nres
4497 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4498 chain_rep(j,chain_length+nres,symetr) &
4499 =chain_rep(j,chain_length+nres,1)
4502 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4504 ! do kkk=1,chain_length
4505 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4509 ! makes copy of chains
4510 write (iout,*) "symetr", symetr
4515 if (symetr.gt.1) then
4522 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4528 write (iout,*) i,icha
4529 do lll=1,chain_length
4531 if (cou.le.nres) then
4533 kupa=mod(lll,chain_length)
4534 iprzes=(kkk-1)*chain_length+lll
4535 if (kupa.eq.0) kupa=chain_length
4536 write (iout,*) "kupa", kupa
4537 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4538 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4545 !-koniec robienia kopii
4548 write (iout,*) "nowa struktura", nperm
4550 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4552 cref(3,i,kkk),cref(1,nres+i,kkk),&
4553 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4555 100 format (//' alpha-carbon coordinates ',&
4556 ' centroid coordinates'/ &
4557 ' ', 6X,'X',11X,'Y',11X,'Z', &
4558 10X,'X',11X,'Y',11X,'Z')
4559 110 format (a,'(',i5,')',6f12.5)
4565 bfrag(i,j)=bfrag(i,j)-ishift
4571 hfrag(i,j)=hfrag(i,j)-ishift
4577 end subroutine readpdb
4578 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4579 !-----------------------------------------------------------------------------
4581 !-----------------------------------------------------------------------------
4582 subroutine read_control
4596 use random, only: random_init
4597 ! implicit real*8 (a-h,o-z)
4598 ! include 'DIMENSIONS'
4600 use prng, only:prng_restart
4602 logical :: OKRandom!, prng_restart
4605 ! include 'COMMON.IOUNITS'
4606 ! include 'COMMON.TIME1'
4607 ! include 'COMMON.THREAD'
4608 ! include 'COMMON.SBRIDGE'
4609 ! include 'COMMON.CONTROL'
4610 ! include 'COMMON.MCM'
4611 ! include 'COMMON.MAP'
4612 ! include 'COMMON.HEADER'
4613 ! include 'COMMON.CSA'
4614 ! include 'COMMON.CHAIN'
4615 ! include 'COMMON.MUCA'
4616 ! include 'COMMON.MD'
4617 ! include 'COMMON.FFIELD'
4618 ! include 'COMMON.INTERACT'
4619 ! include 'COMMON.SETUP'
4620 !el integer :: KDIAG,ICORFL,IXDR
4621 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4622 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4623 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4624 ! character(len=80) :: ucase
4625 character(len=640) :: controlcard
4627 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4633 read (INP,'(a)') titel
4634 call card_concat(controlcard,.true.)
4635 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4636 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4637 call reada(controlcard,'SEED',seed,0.0D0)
4638 call random_init(seed)
4639 ! Set up the time limit (caution! The time must be input in minutes!)
4640 read_cart=index(controlcard,'READ_CART').gt.0
4641 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4642 call readi(controlcard,'SYM',symetr,1)
4643 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4644 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4645 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4646 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4647 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4648 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4649 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4650 call reada(controlcard,'DRMS',drms,0.1D0)
4651 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4652 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4653 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4654 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4655 write (iout,'(a,f10.1)')'DRMS = ',drms
4656 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4657 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4659 call readi(controlcard,'NZ_START',nz_start,0)
4660 call readi(controlcard,'NZ_END',nz_end,0)
4661 call readi(controlcard,'IZ_SC',iz_sc,0)
4662 timlim=60.0D0*timlim
4663 safety = 60.0d0*safety
4666 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4667 !C SHIELD keyword sets if the shielding effect of side-chains is used
4668 !C 0 denots no shielding is used all peptide are equally despite the
4669 !C solvent accesible area
4670 !C 1 the newly introduced function
4671 !C 2 reseved for further possible developement
4672 call readi(controlcard,'SHIELD',shield_mode,0)
4673 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4674 write(iout,*) "shield_mode",shield_mode
4675 call readi(controlcard,'TORMODE',tor_mode,0)
4676 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4677 write(iout,*) "torsional and valence angle mode",tor_mode
4679 !C Varibles set size of box
4680 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4681 protein=index(controlcard,"PROTEIN").gt.0
4682 ions=index(controlcard,"IONS").gt.0
4683 call readi(controlcard,'OLDION',oldion,1)
4684 nucleic=index(controlcard,"NUCLEIC").gt.0
4685 write (iout,*) "with_theta_constr ",with_theta_constr
4686 AFMlog=(index(controlcard,'AFM'))
4687 selfguide=(index(controlcard,'SELFGUIDE'))
4688 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4689 call readi(controlcard,'GENCONSTR',genconstr,0)
4690 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4691 call reada(controlcard,'BOXY',boxysize,100.0d0)
4692 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4693 call readi(controlcard,'TUBEMOD',tubemode,0)
4694 print *,"SCELE",scelemode
4695 call readi(controlcard,"SCELEMODE",scelemode,0)
4696 print *,"SCELE",scelemode
4698 ! elemode = 0 is orignal UNRES electrostatics
4699 ! elemode = 1 is "Momo" potentials in progress
4700 ! elemode = 2 is in development EVALD
4703 write (iout,*) TUBEmode,"TUBEMODE"
4704 if (TUBEmode.gt.0) then
4705 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4706 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4707 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4708 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4709 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4710 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4711 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4712 buftubebot=bordtubebot+tubebufthick
4713 buftubetop=bordtubetop-tubebufthick
4716 ! CUTOFFF ON ELECTROSTATICS
4717 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4718 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4719 write(iout,*) "R_CUT_ELE=",r_cut_ele
4720 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4721 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4722 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4724 ! Lipidec parameters
4725 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4726 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4727 if (lipthick.gt.0.0d0) then
4728 bordliptop=(boxzsize+lipthick)/2.0
4729 bordlipbot=bordliptop-lipthick
4730 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4731 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4732 buflipbot=bordlipbot+lipbufthick
4733 bufliptop=bordliptop-lipbufthick
4734 if ((lipbufthick*2.0d0).gt.lipthick) &
4735 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4736 endif !lipthick.gt.0
4737 write(iout,*) "bordliptop=",bordliptop
4738 write(iout,*) "bordlipbot=",bordlipbot
4739 write(iout,*) "bufliptop=",bufliptop
4740 write(iout,*) "buflipbot=",buflipbot
4741 write (iout,*) "SHIELD MODE",shield_mode
4743 !C-------------------------
4744 minim=(index(controlcard,'MINIMIZE').gt.0)
4745 dccart=(index(controlcard,'CART').gt.0)
4746 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4747 overlapsc=.not.overlapsc
4748 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4749 searchsc=.not.searchsc
4750 sideadd=(index(controlcard,'SIDEADD').gt.0)
4751 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4752 outpdb=(index(controlcard,'PDBOUT').gt.0)
4753 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4754 pdbref=(index(controlcard,'PDBREF').gt.0)
4755 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4756 indpdb=index(controlcard,'PDBSTART')
4757 extconf=(index(controlcard,'EXTCONF').gt.0)
4758 call readi(controlcard,'IPRINT',iprint,0)
4759 call readi(controlcard,'MAXGEN',maxgen,10000)
4760 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4761 call readi(controlcard,"KDIAG",kdiag,0)
4762 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4763 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4764 write (iout,*) "RESCALE_MODE",rescale_mode
4765 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4766 if (index(controlcard,'REGULAR').gt.0.0D0) then
4767 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4771 if (index(controlcard,'CHECKGRAD').gt.0) then
4773 if (index(controlcard,'CART').gt.0) then
4775 elseif (index(controlcard,'CARINT').gt.0) then
4780 elseif (index(controlcard,'THREAD').gt.0) then
4782 call readi(controlcard,'THREAD',nthread,0)
4783 if (nthread.gt.0) then
4784 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4787 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4788 stop 'Error termination in Read_Control.'
4790 else if (index(controlcard,'MCMA').gt.0) then
4792 else if (index(controlcard,'MCEE').gt.0) then
4794 else if (index(controlcard,'MULTCONF').gt.0) then
4796 else if (index(controlcard,'MAP').gt.0) then
4798 call readi(controlcard,'MAP',nmap,0)
4799 else if (index(controlcard,'CSA').gt.0) then
4801 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4803 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4806 !fcm else if (index(controlcard,'MCMF').gt.0) then
4808 else if (index(controlcard,'SOFTREG').gt.0) then
4810 else if (index(controlcard,'CHECK_BOND').gt.0) then
4812 else if (index(controlcard,'TEST').gt.0) then
4814 else if (index(controlcard,'MD').gt.0) then
4816 else if (index(controlcard,'RE ').gt.0) then
4820 lmuca=index(controlcard,'MUCA').gt.0
4821 call readi(controlcard,'MUCADYN',mucadyn,0)
4822 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4823 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4825 write (iout,*) 'MUCADYN=',mucadyn
4826 write (iout,*) 'MUCASMOOTH=',muca_smooth
4829 iscode=index(controlcard,'ONE_LETTER')
4830 indphi=index(controlcard,'PHI')
4831 indback=index(controlcard,'BACK')
4832 iranconf=index(controlcard,'RAND_CONF')
4833 i2ndstr=index(controlcard,'USE_SEC_PRED')
4834 gradout=index(controlcard,'GRADOUT').gt.0
4835 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4836 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4837 if (me.eq.king .or. .not.out1file ) &
4838 write (iout,*) "DISTCHAINMAX",distchainmax
4840 if(me.eq.king.or..not.out1file) &
4841 write (iout,'(2a)') diagmeth(kdiag),&
4842 ' routine used to diagonalize matrices.'
4843 if (shield_mode.gt.0) then
4844 pi=4.0D0*datan(1.0D0)
4845 !C VSolvSphere the volume of solving sphere
4847 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4848 !C there will be no distinction between proline peptide group and normal peptide
4849 !C group in case of shielding parameters
4850 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4851 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4852 write (iout,*) VSolvSphere,VSolvSphere_div
4853 !C long axis of side chain
4855 ! long_r_sidechain(i)=vbldsc0(1,i)
4856 ! short_r_sidechain(i)=sigma0(i)
4857 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4864 end subroutine read_control
4865 !-----------------------------------------------------------------------------
4866 subroutine read_REMDpar
4868 ! Read REMD settings
4875 use control_data, only:out1file
4876 ! implicit real*8 (a-h,o-z)
4877 ! include 'DIMENSIONS'
4878 ! include 'COMMON.IOUNITS'
4879 ! include 'COMMON.TIME1'
4880 ! include 'COMMON.MD'
4883 !el include 'COMMON.LANGEVIN'
4885 !el include 'COMMON.LANGEVIN.lang0'
4887 ! include 'COMMON.INTERACT'
4888 ! include 'COMMON.NAMES'
4889 ! include 'COMMON.GEO'
4890 ! include 'COMMON.REMD'
4891 ! include 'COMMON.CONTROL'
4892 ! include 'COMMON.SETUP'
4893 ! character(len=80) :: ucase
4894 character(len=320) :: controlcard
4895 character(len=3200) :: controlcard1
4896 integer :: iremd_m_total
4899 ! real(kind=8) :: var,ene
4901 if(me.eq.king.or..not.out1file) &
4902 write (iout,*) "REMD setup"
4904 call card_concat(controlcard,.true.)
4905 call readi(controlcard,"NREP",nrep,3)
4906 call readi(controlcard,"NSTEX",nstex,1000)
4907 call reada(controlcard,"RETMIN",retmin,10.0d0)
4908 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4909 mremdsync=(index(controlcard,'SYNC').gt.0)
4910 call readi(controlcard,"NSYN",i_sync_step,100)
4911 restart1file=(index(controlcard,'REST1FILE').gt.0)
4912 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4913 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4914 if(max_cache_traj_use.gt.max_cache_traj) &
4915 max_cache_traj_use=max_cache_traj
4916 if(me.eq.king.or..not.out1file) then
4917 !d if (traj1file) then
4918 !rc caching is in testing - NTWX is not ignored
4919 !d write (iout,*) "NTWX value is ignored"
4920 !d write (iout,*) " trajectory is stored to one file by master"
4921 !d write (iout,*) " before exchange at NSTEX intervals"
4923 write (iout,*) "NREP= ",nrep
4924 write (iout,*) "NSTEX= ",nstex
4925 write (iout,*) "SYNC= ",mremdsync
4926 write (iout,*) "NSYN= ",i_sync_step
4927 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4930 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4931 if (index(controlcard,'TLIST').gt.0) then
4933 call card_concat(controlcard1,.true.)
4934 read(controlcard1,*) (remd_t(i),i=1,nrep)
4935 if(me.eq.king.or..not.out1file) &
4936 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4939 if (index(controlcard,'MLIST').gt.0) then
4941 call card_concat(controlcard1,.true.)
4942 read(controlcard1,*) (remd_m(i),i=1,nrep)
4943 if(me.eq.king.or..not.out1file) then
4944 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4947 iremd_m_total=iremd_m_total+remd_m(i)
4949 write (iout,*) 'Total number of replicas ',iremd_m_total
4952 if(me.eq.king.or..not.out1file) &
4953 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4955 end subroutine read_REMDpar
4956 !-----------------------------------------------------------------------------
4957 subroutine read_MDpar
4961 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
4963 use geometry_data, only: pi
4965 ! implicit real*8 (a-h,o-z)
4966 ! include 'DIMENSIONS'
4967 ! include 'COMMON.IOUNITS'
4968 ! include 'COMMON.TIME1'
4969 ! include 'COMMON.MD'
4972 !el include 'COMMON.LANGEVIN'
4974 !el include 'COMMON.LANGEVIN.lang0'
4976 ! include 'COMMON.INTERACT'
4977 ! include 'COMMON.NAMES'
4978 ! include 'COMMON.GEO'
4979 ! include 'COMMON.SETUP'
4980 ! include 'COMMON.CONTROL'
4981 ! include 'COMMON.SPLITELE'
4982 ! character(len=80) :: ucase
4983 character(len=320) :: controlcard
4988 call card_concat(controlcard,.true.)
4989 call readi(controlcard,"NSTEP",n_timestep,1000000)
4990 call readi(controlcard,"NTWE",ntwe,100)
4991 call readi(controlcard,"NTWX",ntwx,1000)
4992 call reada(controlcard,"DT",d_time,1.0d-1)
4993 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4994 call reada(controlcard,"DAMAX",damax,1.0d1)
4995 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4996 call readi(controlcard,"LANG",lang,0)
4997 RESPA = index(controlcard,"RESPA") .gt. 0
4998 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4999 ntime_split0=ntime_split
5000 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5001 ntime_split0=ntime_split
5002 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5003 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5004 rest = index(controlcard,"REST").gt.0
5005 tbf = index(controlcard,"TBF").gt.0
5006 usampl = index(controlcard,"USAMPL").gt.0
5007 mdpdb = index(controlcard,"MDPDB").gt.0
5008 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5009 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5010 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5011 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5012 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5013 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5014 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5015 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5016 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5017 large = index(controlcard,"LARGE").gt.0
5018 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5019 rattle = index(controlcard,"RATTLE").gt.0
5020 preminim=(index(controlcard,'PREMINIM').gt.0)
5021 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5022 write (iout,*) "PREMINIM ",preminim
5023 dccart=(index(controlcard,'CART').gt.0)
5024 if (preminim) call read_minim
5025 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5031 if(me.eq.king.or..not.out1file) then
5033 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5035 write (iout,'(a)') "The units are:"
5036 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5037 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5038 " acceleration: angstrom/(48.9 fs)**2"
5039 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5041 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5042 write (iout,'(a60,f10.5,a)') &
5043 "Initial time step of numerical integration:",d_time,&
5045 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5047 write (iout,'(2a,i4,a)') &
5048 "A-MTS algorithm used; initial time step for fast-varying",&
5049 " short-range forces split into",ntime_split," steps."
5050 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5051 r_cut," lambda",rlamb
5053 write (iout,'(2a,f10.5)') &
5054 "Maximum acceleration threshold to reduce the time step",&
5055 "/increase split number:",damax
5056 write (iout,'(2a,f10.5)') &
5057 "Maximum predicted energy drift to reduce the timestep",&
5058 "/increase split number:",edriftmax
5059 write (iout,'(a60,f10.5)') &
5060 "Maximum velocity threshold to reduce velocities:",dvmax
5061 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5062 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5063 if (rattle) write (iout,'(a60)') &
5064 "Rattle algorithm used to constrain the virtual bonds"
5068 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5069 call reada(controlcard,"RWAT",rwat,1.4d0)
5070 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5071 surfarea=index(controlcard,"SURFAREA").gt.0
5072 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5073 if(me.eq.king.or..not.out1file)then
5074 write (iout,'(/a,$)') "Langevin dynamics calculation"
5076 write (iout,'(a/)') &
5077 " with direct integration of Langevin equations"
5078 else if (lang.eq.2) then
5079 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5080 else if (lang.eq.3) then
5081 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5082 else if (lang.eq.4) then
5083 write (iout,'(a/)') " in overdamped mode"
5085 write (iout,'(//a,i5)') &
5086 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5089 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5090 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5091 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5092 write (iout,'(a60,f10.5)') &
5093 "Scaling factor of the friction forces:",scal_fric
5094 if (surfarea) write (iout,'(2a,i10,a)') &
5095 "Friction coefficients will be scaled by solvent-accessible",&
5096 " surface area every",reset_fricmat," steps."
5098 ! Calculate friction coefficients and bounds of stochastic forces
5099 eta=6*pi*cPoise*etawat
5100 if(me.eq.king.or..not.out1file) &
5101 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5104 do j=1,5 !types of molecules
5105 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5106 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5108 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5109 do j=1,5 !types of molecules
5111 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5112 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5116 if(me.eq.king.or..not.out1file)then
5118 write (iout,'(/2a/)') &
5119 "Radii of site types and friction coefficients and std's of",&
5120 " stochastic forces of fully exposed sites"
5121 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5124 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5125 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5130 if(me.eq.king.or..not.out1file)then
5131 write (iout,'(a)') "Berendsen bath calculation"
5132 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5133 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5135 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5136 count_reset_moment," steps"
5138 write (iout,'(a,i10,a)') &
5139 "Velocities will be reset at random every",count_reset_vel,&
5143 if(me.eq.king.or..not.out1file) &
5144 write (iout,'(a31)') "Microcanonical mode calculation"
5146 if(me.eq.king.or..not.out1file)then
5147 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5149 write(iout,*) "MD running with constraints."
5150 write(iout,*) "Equilibration time ", eq_time, " mtus."
5151 write(iout,*) "Constraining ", nfrag," fragments."
5152 write(iout,*) "Length of each fragment, weight and q0:"
5154 write (iout,*) "Set of restraints #",iset
5156 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5157 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5159 write(iout,*) "constraints between ", npair, "fragments."
5160 write(iout,*) "constraint pairs, weights and q0:"
5162 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5163 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5165 write(iout,*) "angle constraints within ", nfrag_back,&
5166 "backbone fragments."
5167 write(iout,*) "fragment, weights:"
5169 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5170 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5171 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5174 iset=mod(kolor,nset)+1
5177 if(me.eq.king.or..not.out1file) &
5178 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5180 end subroutine read_MDpar
5181 !-----------------------------------------------------------------------------
5185 ! implicit real*8 (a-h,o-z)
5186 ! include 'DIMENSIONS'
5187 ! include 'COMMON.MAP'
5188 ! include 'COMMON.IOUNITS'
5189 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5190 character(len=80) :: mapcard !,ucase
5193 ! real(kind=8) :: var,ene
5196 read (inp,'(a)') mapcard
5197 mapcard=ucase(mapcard)
5198 if (index(mapcard,'PHI').gt.0) then
5200 else if (index(mapcard,'THE').gt.0) then
5202 else if (index(mapcard,'ALP').gt.0) then
5204 else if (index(mapcard,'OME').gt.0) then
5207 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5208 stop 'Error - illegal variable spec in MAP card.'
5210 call readi (mapcard,'RES1',res1(imap),0)
5211 call readi (mapcard,'RES2',res2(imap),0)
5212 if (res1(imap).eq.0) then
5213 res1(imap)=res2(imap)
5214 else if (res2(imap).eq.0) then
5215 res2(imap)=res1(imap)
5217 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5218 write (iout,'(a)') &
5219 'Error - illegal definition of variable group in MAP.'
5220 stop 'Error - illegal definition of variable group in MAP.'
5222 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5223 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5224 call readi(mapcard,'NSTEP',nstep(imap),0)
5225 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5226 write (iout,'(a)') &
5227 'Illegal boundary and/or step size specification in MAP.'
5228 stop 'Illegal boundary and/or step size specification in MAP.'
5232 end subroutine map_read
5233 !-----------------------------------------------------------------------------
5236 use control_data, only: vdisulf
5238 ! implicit real*8 (a-h,o-z)
5239 ! include 'DIMENSIONS'
5240 ! include 'COMMON.IOUNITS'
5241 ! include 'COMMON.GEO'
5242 ! include 'COMMON.CSA'
5243 ! include 'COMMON.BANK'
5244 ! include 'COMMON.CONTROL'
5245 ! character(len=80) :: ucase
5246 character(len=620) :: mcmcard
5248 ! integer :: ntf,ik,iw_pdb
5249 ! real(kind=8) :: var,ene
5251 call card_concat(mcmcard,.true.)
5253 call readi(mcmcard,'NCONF',nconf,50)
5254 call readi(mcmcard,'NADD',nadd,0)
5255 call readi(mcmcard,'JSTART',jstart,1)
5256 call readi(mcmcard,'JEND',jend,1)
5257 call readi(mcmcard,'NSTMAX',nstmax,500000)
5258 call readi(mcmcard,'N0',n0,1)
5259 call readi(mcmcard,'N1',n1,6)
5260 call readi(mcmcard,'N2',n2,4)
5261 call readi(mcmcard,'N3',n3,0)
5262 call readi(mcmcard,'N4',n4,0)
5263 call readi(mcmcard,'N5',n5,0)
5264 call readi(mcmcard,'N6',n6,10)
5265 call readi(mcmcard,'N7',n7,0)
5266 call readi(mcmcard,'N8',n8,0)
5267 call readi(mcmcard,'N9',n9,0)
5268 call readi(mcmcard,'N14',n14,0)
5269 call readi(mcmcard,'N15',n15,0)
5270 call readi(mcmcard,'N16',n16,0)
5271 call readi(mcmcard,'N17',n17,0)
5272 call readi(mcmcard,'N18',n18,0)
5274 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5276 call readi(mcmcard,'NDIFF',ndiff,2)
5277 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5278 call readi(mcmcard,'IS1',is1,1)
5279 call readi(mcmcard,'IS2',is2,8)
5280 call readi(mcmcard,'NRAN0',nran0,4)
5281 call readi(mcmcard,'NRAN1',nran1,2)
5282 call readi(mcmcard,'IRR',irr,1)
5283 call readi(mcmcard,'NSEED',nseed,20)
5284 call readi(mcmcard,'NTOTAL',ntotal,10000)
5285 call reada(mcmcard,'CUT1',cut1,2.0d0)
5286 call reada(mcmcard,'CUT2',cut2,5.0d0)
5287 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5288 call readi(mcmcard,'ICMAX',icmax,3)
5289 call readi(mcmcard,'IRESTART',irestart,0)
5290 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5293 call reada(mcmcard,'DELE',dele,20.0d0)
5294 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5295 call readi(mcmcard,'IREF',iref,0)
5296 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5297 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5298 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5299 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5300 write (iout,*) "NCONF_IN",nconf_in
5302 end subroutine csaread
5303 !-----------------------------------------------------------------------------
5307 use control_data, only: MaxMoveType
5310 ! implicit real*8 (a-h,o-z)
5311 ! include 'DIMENSIONS'
5312 ! include 'COMMON.MCM'
5313 ! include 'COMMON.MCE'
5314 ! include 'COMMON.IOUNITS'
5315 ! character(len=80) :: ucase
5316 character(len=320) :: mcmcard
5319 ! real(kind=8) :: var,ene
5321 call card_concat(mcmcard,.true.)
5322 call readi(mcmcard,'MAXACC',maxacc,100)
5323 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5324 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5325 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5326 call readi(mcmcard,'MAXREPM',maxrepm,200)
5327 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5328 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5329 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5330 call reada(mcmcard,'E_UP',e_up,5.0D0)
5331 call reada(mcmcard,'DELTE',delte,0.1D0)
5332 call readi(mcmcard,'NSWEEP',nsweep,5)
5333 call readi(mcmcard,'NSTEPH',nsteph,0)
5334 call readi(mcmcard,'NSTEPC',nstepc,0)
5335 call reada(mcmcard,'TMIN',tmin,298.0D0)
5336 call reada(mcmcard,'TMAX',tmax,298.0D0)
5337 call readi(mcmcard,'NWINDOW',nwindow,0)
5338 call readi(mcmcard,'PRINT_MC',print_mc,0)
5339 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5340 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5341 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5342 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5343 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5344 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5345 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5346 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5347 if (nwindow.gt.0) then
5348 allocate(winstart(nwindow)) !!el (maxres)
5349 allocate(winend(nwindow)) !!el
5350 allocate(winlen(nwindow)) !!el
5351 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5353 winlen(i)=winend(i)-winstart(i)+1
5356 if (tmax.lt.tmin) tmax=tmin
5357 if (tmax.eq.tmin) then
5361 if (nstepc.gt.0 .and. nsteph.gt.0) then
5362 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5363 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5365 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5366 ! Probabilities of different move types
5367 sumpro_type(0)=0.0D0
5368 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5369 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5370 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5371 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5372 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5373 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5374 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5376 print *,'i',i,' sumprotype',sumpro_type(i)
5377 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5378 print *,'i',i,' sumprotype',sumpro_type(i)
5381 end subroutine mcmread
5382 !-----------------------------------------------------------------------------
5383 subroutine read_minim
5386 ! implicit real*8 (a-h,o-z)
5387 ! include 'DIMENSIONS'
5388 ! include 'COMMON.MINIM'
5389 ! include 'COMMON.IOUNITS'
5390 ! character(len=80) :: ucase
5391 character(len=320) :: minimcard
5393 ! integer :: ntf,ik,iw_pdb
5394 ! real(kind=8) :: var,ene
5396 call card_concat(minimcard,.true.)
5397 call readi(minimcard,'MAXMIN',maxmin,2000)
5398 call readi(minimcard,'MAXFUN',maxfun,5000)
5399 call readi(minimcard,'MINMIN',minmin,maxmin)
5400 call readi(minimcard,'MINFUN',minfun,maxmin)
5401 call reada(minimcard,'TOLF',tolf,1.0D-2)
5402 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5403 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5404 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5405 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5406 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5407 'Options in energy minimization:'
5408 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5409 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5410 'MinMin:',MinMin,' MinFun:',MinFun,&
5411 ' TolF:',TolF,' RTolF:',RTolF
5413 end subroutine read_minim
5414 !-----------------------------------------------------------------------------
5415 subroutine openunits
5417 use MD_data, only: usampl
5420 use control_data, only:out1file
5421 use control, only: getenv_loc
5422 ! implicit real*8 (a-h,o-z)
5423 ! include 'DIMENSIONS'
5426 character(len=16) :: form,nodename
5427 integer :: nodelen,ierror,npos
5429 ! include 'COMMON.SETUP'
5430 ! include 'COMMON.IOUNITS'
5431 ! include 'COMMON.MD'
5432 ! include 'COMMON.CONTROL'
5433 integer :: lenpre,lenpot,lentmp !,ilen
5435 character(len=3) :: out1file_text !,ucase
5436 character(len=3) :: ll
5439 ! integer :: ntf,ik,iw_pdb
5440 ! real(kind=8) :: var,ene
5442 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5443 call getenv_loc("PREFIX",prefix)
5445 call getenv_loc("POT",pot)
5446 call getenv_loc("DIRTMP",tmpdir)
5447 call getenv_loc("CURDIR",curdir)
5448 call getenv_loc("OUT1FILE",out1file_text)
5449 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5450 out1file_text=ucase(out1file_text)
5451 if (out1file_text(1:1).eq."Y") then
5454 out1file=fg_rank.gt.0
5459 if (lentmp.gt.0) then
5460 write (*,'(80(1h!))')
5461 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5462 write (*,'(80(1h!))')
5463 write (*,*)"All output files will be on node /tmp directory."
5465 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5466 if (me.eq.king) then
5467 write (*,*) "The master node is ",nodename
5468 else if (fg_rank.eq.0) then
5469 write (*,*) "I am the CG slave node ",nodename
5471 write (*,*) "I am the FG slave node ",nodename
5474 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5475 lenpre = lentmp+lenpre+1
5477 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5478 ! Get the names and open the input files
5479 #if defined(WINIFL) || defined(WINPGI)
5480 open(1,file=pref_orig(:ilen(pref_orig))// &
5481 '.inp',status='old',readonly,shared)
5482 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5483 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5484 ! Get parameter filenames and open the parameter files.
5485 call getenv_loc('BONDPAR',bondname)
5486 open (ibond,file=bondname,status='old',readonly,shared)
5487 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5488 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5489 call getenv_loc('THETPAR',thetname)
5490 open (ithep,file=thetname,status='old',readonly,shared)
5491 call getenv_loc('ROTPAR',rotname)
5492 open (irotam,file=rotname,status='old',readonly,shared)
5493 call getenv_loc('TORPAR',torname)
5494 open (itorp,file=torname,status='old',readonly,shared)
5495 call getenv_loc('TORDPAR',tordname)
5496 open (itordp,file=tordname,status='old',readonly,shared)
5497 call getenv_loc('FOURIER',fouriername)
5498 open (ifourier,file=fouriername,status='old',readonly,shared)
5499 call getenv_loc('ELEPAR',elename)
5500 open (ielep,file=elename,status='old',readonly,shared)
5501 call getenv_loc('SIDEPAR',sidename)
5502 open (isidep,file=sidename,status='old',readonly,shared)
5504 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5505 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5506 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5507 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5508 call getenv_loc('TORPAR_NUCL',torname_nucl)
5509 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5510 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5511 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5512 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5513 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5516 #elif (defined CRAY) || (defined AIX)
5517 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5519 ! print *,"Processor",myrank," opened file 1"
5520 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5521 ! print *,"Processor",myrank," opened file 9"
5522 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5523 ! Get parameter filenames and open the parameter files.
5524 call getenv_loc('BONDPAR',bondname)
5525 open (ibond,file=bondname,status='old',action='read')
5526 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5527 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5529 ! print *,"Processor",myrank," opened file IBOND"
5530 call getenv_loc('THETPAR',thetname)
5531 open (ithep,file=thetname,status='old',action='read')
5532 ! print *,"Processor",myrank," opened file ITHEP"
5533 call getenv_loc('ROTPAR',rotname)
5534 open (irotam,file=rotname,status='old',action='read')
5535 ! print *,"Processor",myrank," opened file IROTAM"
5536 call getenv_loc('TORPAR',torname)
5537 open (itorp,file=torname,status='old',action='read')
5538 ! print *,"Processor",myrank," opened file ITORP"
5539 call getenv_loc('TORDPAR',tordname)
5540 open (itordp,file=tordname,status='old',action='read')
5541 ! print *,"Processor",myrank," opened file ITORDP"
5542 call getenv_loc('SCCORPAR',sccorname)
5543 open (isccor,file=sccorname,status='old',action='read')
5544 ! print *,"Processor",myrank," opened file ISCCOR"
5545 call getenv_loc('FOURIER',fouriername)
5546 open (ifourier,file=fouriername,status='old',action='read')
5547 ! print *,"Processor",myrank," opened file IFOURIER"
5548 call getenv_loc('ELEPAR',elename)
5549 open (ielep,file=elename,status='old',action='read')
5550 ! print *,"Processor",myrank," opened file IELEP"
5551 call getenv_loc('SIDEPAR',sidename)
5552 open (isidep,file=sidename,status='old',action='read')
5554 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5555 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5556 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5557 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5558 call getenv_loc('TORPAR_NUCL',torname_nucl)
5559 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5560 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5561 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5562 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5563 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5565 call getenv_loc('LIPTRANPAR',liptranname)
5566 open (iliptranpar,file=liptranname,status='old',action='read')
5567 call getenv_loc('TUBEPAR',tubename)
5568 open (itube,file=tubename,status='old',action='read')
5569 call getenv_loc('IONPAR',ionname)
5570 open (iion,file=ionname,status='old',action='read')
5571 call getenv_loc('IONPAR_TRAN',iontranname)
5572 open (iiontran,file=iontranname,status='old',action='read')
5573 ! print *,"Processor",myrank," opened file ISIDEP"
5574 ! print *,"Processor",myrank," opened parameter files"
5576 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5577 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5578 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5579 ! Get parameter filenames and open the parameter files.
5580 call getenv_loc('BONDPAR',bondname)
5581 open (ibond,file=bondname,status='old')
5582 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5583 open (ibond_nucl,file=bondname_nucl,status='old')
5585 call getenv_loc('THETPAR',thetname)
5586 open (ithep,file=thetname,status='old')
5587 call getenv_loc('ROTPAR',rotname)
5588 open (irotam,file=rotname,status='old')
5589 call getenv_loc('TORPAR',torname)
5590 open (itorp,file=torname,status='old')
5591 call getenv_loc('TORDPAR',tordname)
5592 open (itordp,file=tordname,status='old')
5593 call getenv_loc('SCCORPAR',sccorname)
5594 open (isccor,file=sccorname,status='old')
5595 call getenv_loc('FOURIER',fouriername)
5596 open (ifourier,file=fouriername,status='old')
5597 call getenv_loc('ELEPAR',elename)
5598 open (ielep,file=elename,status='old')
5599 call getenv_loc('SIDEPAR',sidename)
5600 open (isidep,file=sidename,status='old')
5602 open (ithep_nucl,file=thetname_nucl,status='old')
5603 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5604 open (irotam_nucl,file=rotname_nucl,status='old')
5605 call getenv_loc('TORPAR_NUCL',torname_nucl)
5606 open (itorp_nucl,file=torname_nucl,status='old')
5607 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5608 open (itordp_nucl,file=tordname_nucl,status='old')
5609 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5610 open (isidep_nucl,file=sidename_nucl,status='old')
5612 call getenv_loc('LIPTRANPAR',liptranname)
5613 open (iliptranpar,file=liptranname,status='old')
5614 call getenv_loc('TUBEPAR',tubename)
5615 open (itube,file=tubename,status='old')
5616 call getenv_loc('IONPAR',ionname)
5617 open (iion,file=ionname,status='old')
5618 call getenv_loc('IONPAR_NUCL',ionnuclname)
5619 open (iionnucl,file=ionnuclname,status='old')
5620 call getenv_loc('IONPAR_TRAN',iontranname)
5621 open (iiontran,file=iontranname,status='old')
5623 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5625 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5626 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5627 ! Get parameter filenames and open the parameter files.
5628 call getenv_loc('BONDPAR',bondname)
5629 open (ibond,file=bondname,status='old',action='read')
5630 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5631 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5632 call getenv_loc('THETPAR',thetname)
5633 open (ithep,file=thetname,status='old',action='read')
5634 call getenv_loc('ROTPAR',rotname)
5635 open (irotam,file=rotname,status='old',action='read')
5636 call getenv_loc('TORPAR',torname)
5637 open (itorp,file=torname,status='old',action='read')
5638 call getenv_loc('TORDPAR',tordname)
5639 open (itordp,file=tordname,status='old',action='read')
5640 call getenv_loc('SCCORPAR',sccorname)
5641 open (isccor,file=sccorname,status='old',action='read')
5643 call getenv_loc('THETPARPDB',thetname_pdb)
5644 print *,"thetname_pdb ",thetname_pdb
5645 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5646 print *,ithep_pdb," opened"
5648 call getenv_loc('FOURIER',fouriername)
5649 open (ifourier,file=fouriername,status='old',readonly)
5650 call getenv_loc('ELEPAR',elename)
5651 open (ielep,file=elename,status='old',readonly)
5652 call getenv_loc('SIDEPAR',sidename)
5653 open (isidep,file=sidename,status='old',readonly)
5655 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5656 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5657 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5658 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5659 call getenv_loc('TORPAR_NUCL',torname_nucl)
5660 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5661 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5662 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5663 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5664 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5665 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5666 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5667 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5668 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5669 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5670 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5671 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5672 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5675 call getenv_loc('LIPTRANPAR',liptranname)
5676 open (iliptranpar,file=liptranname,status='old',action='read')
5677 call getenv_loc('LIPBOND',lipbondname)
5678 open (ilipbond,file=lipbondname,status='old',action='read')
5679 call getenv_loc('LIPNONBOND',lipnonbondname)
5680 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5681 call getenv_loc('TUBEPAR',tubename)
5682 open (itube,file=tubename,status='old',action='read')
5683 call getenv_loc('IONPAR',ionname)
5684 open (iion,file=ionname,status='old',action='read')
5685 call getenv_loc('IONPAR_NUCL',ionnuclname)
5686 open (iionnucl,file=ionnuclname,status='old',action='read')
5687 call getenv_loc('IONPAR_TRAN',iontranname)
5688 open (iiontran,file=iontranname,status='old',action='read')
5692 call getenv_loc('ROTPARPDB',rotname_pdb)
5693 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5696 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5697 #if defined(WINIFL) || defined(WINPGI)
5698 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5699 #elif (defined CRAY) || (defined AIX)
5700 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5702 open (iscpp_nucl,file=scpname_nucl,status='old')
5704 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5709 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5710 ! Use -DOLDSCP to use hard-coded constants instead.
5712 call getenv_loc('SCPPAR',scpname)
5713 #if defined(WINIFL) || defined(WINPGI)
5714 open (iscpp,file=scpname,status='old',readonly,shared)
5715 #elif (defined CRAY) || (defined AIX)
5716 open (iscpp,file=scpname,status='old',action='read')
5718 open (iscpp,file=scpname,status='old')
5720 open (iscpp,file=scpname,status='old',action='read')
5723 call getenv_loc('PATTERN',patname)
5724 #if defined(WINIFL) || defined(WINPGI)
5725 open (icbase,file=patname,status='old',readonly,shared)
5726 #elif (defined CRAY) || (defined AIX)
5727 open (icbase,file=patname,status='old',action='read')
5729 open (icbase,file=patname,status='old')
5731 open (icbase,file=patname,status='old',action='read')
5734 ! Open output file only for CG processes
5735 ! print *,"Processor",myrank," fg_rank",fg_rank
5736 if (fg_rank.eq.0) then
5738 if (nodes.eq.1) then
5741 npos = dlog10(dfloat(nodes-1))+1
5743 if (npos.lt.3) npos=3
5744 write (liczba,'(i1)') npos
5745 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5747 write (liczba,form) me
5748 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5749 liczba(:ilen(liczba))
5750 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5752 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5754 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5755 liczba(:ilen(liczba))//'.mol2'
5756 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5757 liczba(:ilen(liczba))//'.stat'
5759 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5760 //liczba(:ilen(liczba))//'.stat')
5761 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5764 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5765 liczba(:ilen(liczba))//'.const'
5770 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5771 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5772 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5773 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5774 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5776 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5778 rest2name=prefix(:ilen(prefix))//'.rst'
5780 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5783 #if defined(AIX) || defined(PGI)
5784 if (me.eq.king .or. .not. out1file) &
5785 open(iout,file=outname,status='unknown')
5787 if (fg_rank.gt.0) then
5788 write (liczba,'(i3.3)') myrank/nfgtasks
5789 write (ll,'(bz,i3.3)') fg_rank
5790 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5795 open(igeom,file=intname,status='unknown',position='append')
5796 open(ipdb,file=pdbname,status='unknown')
5797 open(imol2,file=mol2name,status='unknown')
5798 open(istat,file=statname,status='unknown',position='append')
5800 !1out open(iout,file=outname,status='unknown')
5803 if (me.eq.king .or. .not.out1file) &
5804 open(iout,file=outname,status='unknown')
5806 if (fg_rank.gt.0) then
5807 write (liczba,'(i3.3)') myrank/nfgtasks
5808 write (ll,'(bz,i3.3)') fg_rank
5809 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5814 open(igeom,file=intname,status='unknown',access='append')
5815 open(ipdb,file=pdbname,status='unknown')
5816 open(imol2,file=mol2name,status='unknown')
5817 open(istat,file=statname,status='unknown',access='append')
5819 !1out open(iout,file=outname,status='unknown')
5822 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5823 csa_seed=prefix(:lenpre)//'.CSA.seed'
5824 csa_history=prefix(:lenpre)//'.CSA.history'
5825 csa_bank=prefix(:lenpre)//'.CSA.bank'
5826 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5827 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5828 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5829 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5830 csa_int=prefix(:lenpre)//'.int'
5831 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5832 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5833 csa_in=prefix(:lenpre)//'.CSA.in'
5834 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5837 write (iout,'(80(1h-))')
5838 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5839 write (iout,'(80(1h-))')
5840 write (iout,*) "Input file : ",&
5841 pref_orig(:ilen(pref_orig))//'.inp'
5842 write (iout,*) "Output file : ",&
5843 outname(:ilen(outname))
5845 write (iout,*) "Sidechain potential file : ",&
5846 sidename(:ilen(sidename))
5848 write (iout,*) "SCp potential file : ",&
5849 scpname(:ilen(scpname))
5851 write (iout,*) "Electrostatic potential file : ",&
5852 elename(:ilen(elename))
5853 write (iout,*) "Cumulant coefficient file : ",&
5854 fouriername(:ilen(fouriername))
5855 write (iout,*) "Torsional parameter file : ",&
5856 torname(:ilen(torname))
5857 write (iout,*) "Double torsional parameter file : ",&
5858 tordname(:ilen(tordname))
5859 write (iout,*) "SCCOR parameter file : ",&
5860 sccorname(:ilen(sccorname))
5861 write (iout,*) "Bond & inertia constant file : ",&
5862 bondname(:ilen(bondname))
5863 write (iout,*) "Bending parameter file : ",&
5864 thetname(:ilen(thetname))
5865 write (iout,*) "Rotamer parameter file : ",&
5866 rotname(:ilen(rotname))
5869 write (iout,*) "Thetpdb parameter file : ",&
5870 thetname_pdb(:ilen(thetname_pdb))
5873 write (iout,*) "Threading database : ",&
5874 patname(:ilen(patname))
5876 write (iout,*)" DIRTMP : ",&
5878 write (iout,'(80(1h-))')
5881 end subroutine openunits
5882 !-----------------------------------------------------------------------------
5885 use geometry_data, only: nres,dc
5887 ! implicit real*8 (a-h,o-z)
5888 ! include 'DIMENSIONS'
5889 ! include 'COMMON.CHAIN'
5890 ! include 'COMMON.IOUNITS'
5891 ! include 'COMMON.MD'
5894 ! real(kind=8) :: var,ene
5896 open(irest2,file=rest2name,status='unknown')
5897 read(irest2,*) totT,EK,potE,totE,t_bath
5900 ! AL 4/17/17: Now reading d_t(0,:) too
5902 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5905 ! AL 4/17/17: Now reading d_c(0,:) too
5907 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5910 read (irest2,*) iset
5914 end subroutine readrst
5915 !-----------------------------------------------------------------------------
5916 subroutine read_fragments
5920 use control_data, only:out1file
5923 ! implicit real*8 (a-h,o-z)
5924 ! include 'DIMENSIONS'
5928 ! include 'COMMON.SETUP'
5929 ! include 'COMMON.CHAIN'
5930 ! include 'COMMON.IOUNITS'
5931 ! include 'COMMON.MD'
5932 ! include 'COMMON.CONTROL'
5935 ! real(kind=8) :: var,ene
5937 read(inp,*) nset,nfrag,npair,nfrag_back
5939 !el from module energy
5940 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5941 if(.not.allocated(wfrag_back)) then
5942 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5943 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5945 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5946 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5948 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5949 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5952 if(me.eq.king.or..not.out1file) &
5953 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5954 " nfrag_back",nfrag_back
5956 read(inp,*) mset(iset)
5958 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5960 if(me.eq.king.or..not.out1file) &
5961 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5962 ifrag(2,i,iset), qinfrag(i,iset)
5965 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5967 if(me.eq.king.or..not.out1file) &
5968 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5969 ipair(2,i,iset), qinpair(i,iset)
5972 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5973 wfrag_back(3,i,iset),&
5974 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5975 if(me.eq.king.or..not.out1file) &
5976 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5977 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5981 end subroutine read_fragments
5982 !-----------------------------------------------------------------------------
5984 !-----------------------------------------------------------------------------
5988 ! implicit real*8 (a-h,o-z)
5989 ! include 'DIMENSIONS'
5990 ! include 'COMMON.CSA'
5991 ! include 'COMMON.BANK'
5992 ! include 'COMMON.IOUNITS'
5994 ! integer :: ntf,ik,iw_pdb
5995 ! real(kind=8) :: var,ene
5997 open(icsa_in,file=csa_in,status="old",err=100)
5998 read(icsa_in,*) nconf
5999 read(icsa_in,*) jstart,jend
6000 read(icsa_in,*) nstmax
6001 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6002 read(icsa_in,*) nran0,nran1,irr
6003 read(icsa_in,*) nseed
6004 read(icsa_in,*) ntotal,cut1,cut2
6005 read(icsa_in,*) estop
6006 read(icsa_in,*) icmax,irestart
6007 read(icsa_in,*) ntbankm,dele,difcut
6008 read(icsa_in,*) iref,rmscut,pnccut
6009 read(icsa_in,*) ndiff
6016 end subroutine csa_read
6017 !-----------------------------------------------------------------------------
6018 subroutine initial_write
6021 ! implicit real*8 (a-h,o-z)
6022 ! include 'DIMENSIONS'
6023 ! include 'COMMON.CSA'
6024 ! include 'COMMON.BANK'
6025 ! include 'COMMON.IOUNITS'
6027 ! integer :: ntf,ik,iw_pdb
6028 ! real(kind=8) :: var,ene
6030 open(icsa_seed,file=csa_seed,status="unknown")
6031 write(icsa_seed,*) "seed"
6033 #if defined(AIX) || defined(PGI)
6034 open(icsa_history,file=csa_history,status="unknown",&
6037 open(icsa_history,file=csa_history,status="unknown",&
6040 write(icsa_history,*) nconf
6041 write(icsa_history,*) jstart,jend
6042 write(icsa_history,*) nstmax
6043 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6044 write(icsa_history,*) nran0,nran1,irr
6045 write(icsa_history,*) nseed
6046 write(icsa_history,*) ntotal,cut1,cut2
6047 write(icsa_history,*) estop
6048 write(icsa_history,*) icmax,irestart
6049 write(icsa_history,*) ntbankm,dele,difcut
6050 write(icsa_history,*) iref,rmscut,pnccut
6051 write(icsa_history,*) ndiff
6053 write(icsa_history,*)
6056 open(icsa_bank1,file=csa_bank1,status="unknown")
6057 write(icsa_bank1,*) 0
6061 end subroutine initial_write
6062 !-----------------------------------------------------------------------------
6063 subroutine restart_write
6066 ! implicit real*8 (a-h,o-z)
6067 ! include 'DIMENSIONS'
6068 ! include 'COMMON.IOUNITS'
6069 ! include 'COMMON.CSA'
6070 ! include 'COMMON.BANK'
6072 ! integer :: ntf,ik,iw_pdb
6073 ! real(kind=8) :: var,ene
6075 #if defined(AIX) || defined(PGI)
6076 open(icsa_history,file=csa_history,position="append")
6078 open(icsa_history,file=csa_history,access="append")
6080 write(icsa_history,*)
6081 write(icsa_history,*) "This is restart"
6082 write(icsa_history,*)
6083 write(icsa_history,*) nconf
6084 write(icsa_history,*) jstart,jend
6085 write(icsa_history,*) nstmax
6086 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6087 write(icsa_history,*) nran0,nran1,irr
6088 write(icsa_history,*) nseed
6089 write(icsa_history,*) ntotal,cut1,cut2
6090 write(icsa_history,*) estop
6091 write(icsa_history,*) icmax,irestart
6092 write(icsa_history,*) ntbankm,dele,difcut
6093 write(icsa_history,*) iref,rmscut,pnccut
6094 write(icsa_history,*) ndiff
6095 write(icsa_history,*)
6096 write(icsa_history,*) "irestart is: ", irestart
6098 write(icsa_history,*)
6102 end subroutine restart_write
6103 !-----------------------------------------------------------------------------
6105 !-----------------------------------------------------------------------------
6106 subroutine write_pdb(npdb,titelloc,ee)
6108 ! implicit real*8 (a-h,o-z)
6109 ! include 'DIMENSIONS'
6110 ! include 'COMMON.IOUNITS'
6111 character(len=50) :: titelloc1
6112 character*(*) :: titelloc
6113 character(len=3) :: zahl
6114 character(len=5) :: liczba5
6116 integer :: npdb !,ilen
6120 ! real(kind=8) :: var,ene
6124 if (npdb.lt.1000) then
6125 call numstr(npdb,zahl)
6126 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6128 if (npdb.lt.10000) then
6129 write(liczba5,'(i1,i4)') 0,npdb
6131 write(liczba5,'(i5)') npdb
6133 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6135 call pdbout(ee,titelloc1,ipdb)
6138 end subroutine write_pdb
6139 !-----------------------------------------------------------------------------
6141 !-----------------------------------------------------------------------------
6142 subroutine write_thread_summary
6143 ! Thread the sequence through a database of known structures
6144 use control_data, only: refstr
6146 use energy_data, only: n_ene_comp
6148 ! implicit real*8 (a-h,o-z)
6149 ! include 'DIMENSIONS'
6151 use MPI_data !include 'COMMON.INFO'
6154 ! include 'COMMON.CONTROL'
6155 ! include 'COMMON.CHAIN'
6156 ! include 'COMMON.DBASE'
6157 ! include 'COMMON.INTERACT'
6158 ! include 'COMMON.VAR'
6159 ! include 'COMMON.THREAD'
6160 ! include 'COMMON.FFIELD'
6161 ! include 'COMMON.SBRIDGE'
6162 ! include 'COMMON.HEADER'
6163 ! include 'COMMON.NAMES'
6164 ! include 'COMMON.IOUNITS'
6165 ! include 'COMMON.TIME1'
6167 integer,dimension(maxthread) :: ip
6168 real(kind=8),dimension(0:n_ene) :: energia
6170 integer :: i,j,ii,jj,ipj,ik,kk,ist
6171 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6173 write (iout,'(30x,a/)') &
6174 ' *********** Summary threading statistics ************'
6175 write (iout,'(a)') 'Initial energies:'
6176 write (iout,'(a4,2x,a12,14a14,3a8)') &
6177 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6178 'RMSnat','NatCONT','NNCONT','RMS'
6179 ! Energy sort patterns
6184 enet=ener(n_ene-1,ip(i))
6187 if (ener(n_ene-1,ip(j)).lt.enet) then
6189 enet=ener(n_ene-1,ip(j))
6201 ist=nres_base(2,ii)+ipatt(2,i)
6203 energia(i)=ener0(kk,i)
6205 etot=ener0(n_ene_comp+1,i)
6206 rmsnat=ener0(n_ene_comp+2,i)
6207 rms=ener0(n_ene_comp+3,i)
6208 frac=ener0(n_ene_comp+4,i)
6209 frac_nn=ener0(n_ene_comp+5,i)
6212 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6213 i,str_nam(ii),ist+1,&
6214 (energia(print_order(kk)),kk=1,nprint_ene),&
6215 etot,rmsnat,frac,frac_nn,rms
6217 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6218 i,str_nam(ii),ist+1,&
6219 (energia(print_order(kk)),kk=1,nprint_ene),etot
6222 write (iout,'(//a)') 'Final energies:'
6223 write (iout,'(a4,2x,a12,17a14,3a8)') &
6224 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6225 'RMSnat','NatCONT','NNCONT','RMS'
6229 ist=nres_base(2,ii)+ipatt(2,i)
6231 energia(kk)=ener(kk,ik)
6233 etot=ener(n_ene_comp+1,i)
6234 rmsnat=ener(n_ene_comp+2,i)
6235 rms=ener(n_ene_comp+3,i)
6236 frac=ener(n_ene_comp+4,i)
6237 frac_nn=ener(n_ene_comp+5,i)
6238 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6239 i,str_nam(ii),ist+1,&
6240 (energia(print_order(kk)),kk=1,nprint_ene),&
6241 etot,rmsnat,frac,frac_nn,rms
6243 write (iout,'(/a/)') 'IEXAM array:'
6244 write (iout,'(i5)') nexcl
6246 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6248 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6249 'Max. time for threading step ',max_time_for_thread,&
6250 'Average time for threading step: ',ave_time_for_thread
6252 end subroutine write_thread_summary
6253 !-----------------------------------------------------------------------------
6254 subroutine write_stat_thread(ithread,ipattern,ist)
6256 use energy_data, only: n_ene_comp
6258 ! implicit real*8 (a-h,o-z)
6259 ! include "DIMENSIONS"
6260 ! include "COMMON.CONTROL"
6261 ! include "COMMON.IOUNITS"
6262 ! include "COMMON.THREAD"
6263 ! include "COMMON.FFIELD"
6264 ! include "COMMON.DBASE"
6265 ! include "COMMON.NAMES"
6266 real(kind=8),dimension(0:n_ene) :: energia
6268 integer :: ithread,ipattern,ist,i
6269 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6271 #if defined(AIX) || defined(PGI)
6272 open(istat,file=statname,position='append')
6274 open(istat,file=statname,access='append')
6277 energia(i)=ener(i,ithread)
6279 etot=ener(n_ene_comp+1,ithread)
6280 rmsnat=ener(n_ene_comp+2,ithread)
6281 rms=ener(n_ene_comp+3,ithread)
6282 frac=ener(n_ene_comp+4,ithread)
6283 frac_nn=ener(n_ene_comp+5,ithread)
6284 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6285 ithread,str_nam(ipattern),ist+1,&
6286 (energia(print_order(i)),i=1,nprint_ene),&
6287 etot,rmsnat,frac,frac_nn,rms
6290 end subroutine write_stat_thread
6291 !-----------------------------------------------------------------------------
6293 !-----------------------------------------------------------------------------
6294 end module io_config