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
534 write (iout,*) 'FTORS factor =',ftors
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
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: buse
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii,ncatprotparm
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
780 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
781 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(nbondterm(ntyp)) !(ntyp)
783 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
784 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(msc(ntyp+1,5)) !(ntyp+1)
786 allocate(isc(ntyp+1,5)) !(ntyp+1)
787 allocate(restok(ntyp+1,5)) !(ntyp+1)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 read (ibond,*) vbldp0,akp,mp,ip,pstok
798 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
799 dsc(i) = vbldsc0(1,i)
803 dsc_inv(i)=1.0D0/dsc(i)
807 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
809 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
810 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
811 dsc(i) = vbldsc0(1,i)
815 dsc_inv(i)=1.0D0/dsc(i)
819 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
822 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
823 ! dsc(i) = vbldsc0_nucl(1,i)
827 ! dsc_inv(i)=1.0D0/dsc(i)
830 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
831 ! do i=1,ntyp_molec(2)
832 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
833 ! aksc_nucl(j,i),abond0_nucl(j,i),&
834 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
835 ! dsc(i) = vbldsc0(1,i)
839 ! dsc_inv(i)=1.0D0/dsc(i)
844 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
845 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
847 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
849 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
850 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
852 write (iout,'(13x,3f10.5)') &
853 vbldsc0(j,i),aksc(j,i),abond0(j,i)
858 read(iion,*) msc(i,5),restok(i,5)
859 print *,msc(i,5),restok(i,5)
863 read (iion,*) ncatprotparm
864 allocate(catprm(ncatprotparm,4))
866 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
869 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
870 !----------------------------------------------------
871 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
872 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
873 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
874 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
875 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
876 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
886 allocate(liptranene(ntyp))
887 !C reading lipid parameters
888 write (iout,*) "iliptranpar",iliptranpar
890 read(iliptranpar,*) pepliptran
893 read(iliptranpar,*) liptranene(i)
894 print *,liptranene(i)
900 ! Read the parameters of the probability distribution/energy expression
901 ! of the virtual-bond valence angles theta
904 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
905 (bthet(j,i,1,1),j=1,2)
906 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
907 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
908 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
912 athet(1,i,1,-1)=athet(1,i,1,1)
913 athet(2,i,1,-1)=athet(2,i,1,1)
914 bthet(1,i,1,-1)=-bthet(1,i,1,1)
915 bthet(2,i,1,-1)=-bthet(2,i,1,1)
916 athet(1,i,-1,1)=-athet(1,i,1,1)
917 athet(2,i,-1,1)=-athet(2,i,1,1)
918 bthet(1,i,-1,1)=bthet(1,i,1,1)
919 bthet(2,i,-1,1)=bthet(2,i,1,1)
923 athet(1,i,-1,-1)=athet(1,-i,1,1)
924 athet(2,i,-1,-1)=-athet(2,-i,1,1)
925 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
926 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
927 athet(1,i,-1,1)=athet(1,-i,1,1)
928 athet(2,i,-1,1)=-athet(2,-i,1,1)
929 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
930 bthet(2,i,-1,1)=bthet(2,-i,1,1)
931 athet(1,i,1,-1)=-athet(1,-i,1,1)
932 athet(2,i,1,-1)=athet(2,-i,1,1)
933 bthet(1,i,1,-1)=bthet(1,-i,1,1)
934 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
939 polthet(j,i)=polthet(j,-i)
942 gthet(j,i)=gthet(j,-i)
950 'Parameters of the virtual-bond valence angles:'
951 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
952 ' ATHETA0 ',' A1 ',' A2 ',&
955 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
956 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
958 write (iout,'(/a/9x,5a/79(1h-))') &
959 'Parameters of the expression for sigma(theta_c):',&
960 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
961 ' ALPH3 ',' SIGMA0C '
963 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
964 (polthet(j,i),j=0,3),sigc0(i)
966 write (iout,'(/a/9x,5a/79(1h-))') &
967 'Parameters of the second gaussian:',&
968 ' THETA0 ',' SIGMA0 ',' G1 ',&
971 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
972 sig0(i),(gthet(j,i),j=1,3)
976 'Parameters of the virtual-bond valence angles:'
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Coefficients of expansion',&
979 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
980 ' b1*10^1 ',' b2*10^1 '
982 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
983 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
984 (10*bthet(j,i,1,1),j=1,2)
986 write (iout,'(/a/9x,5a/79(1h-))') &
987 'Parameters of the expression for sigma(theta_c):',&
988 ' alpha0 ',' alph1 ',' alph2 ',&
989 ' alhp3 ',' sigma0c '
991 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
992 (polthet(j,i),j=0,3),sigc0(i)
994 write (iout,'(/a/9x,5a/79(1h-))') &
995 'Parameters of the second gaussian:',&
996 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
999 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1000 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1006 ! Read the parameters of Utheta determined from ab initio surfaces
1007 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1009 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1010 ntheterm3,nsingle,ndouble
1011 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1013 !----------------------------------------------------
1014 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1015 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1016 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1017 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1018 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1019 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1020 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1021 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1022 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1023 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1024 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1025 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1026 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1029 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1030 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1031 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1035 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1037 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1039 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1041 ithetyp(i)=-ithetyp(-i)
1044 aa0thet(:,:,:,:)=0.0d0
1045 aathet(:,:,:,:,:)=0.0d0
1046 bbthet(:,:,:,:,:,:)=0.0d0
1047 ccthet(:,:,:,:,:,:)=0.0d0
1048 ddthet(:,:,:,:,:,:)=0.0d0
1049 eethet(:,:,:,:,:,:)=0.0d0
1050 ffthet(:,:,:,:,:,:,:)=0.0d0
1051 ggthet(:,:,:,:,:,:,:)=0.0d0
1053 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1055 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1056 ! VAR:1=non-glicyne non-proline 2=proline
1057 ! VAR:negative values for D-aminoacid
1059 do j=-nthetyp,nthetyp
1060 do k=-nthetyp,nthetyp
1061 read (ithep,'(6a)',end=111,err=111) res1
1062 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1063 ! VAR: aa0thet is variable describing the average value of Foureir
1064 ! VAR: expansion series
1065 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1066 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1067 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1068 read (ithep,*,end=111,err=111) &
1069 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1070 read (ithep,*,end=111,err=111) &
1071 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1072 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1073 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1074 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1076 read (ithep,*,end=111,err=111) &
1077 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1078 ffthet(lll,llll,ll,i,j,k,iblock),&
1079 ggthet(llll,lll,ll,i,j,k,iblock),&
1080 ggthet(lll,llll,ll,i,j,k,iblock),&
1081 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1086 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1087 ! coefficients of theta-and-gamma-dependent terms are zero.
1088 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1089 ! RECOMENTDED AFTER VERSION 3.3)
1093 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1094 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1096 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1097 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1100 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1102 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1105 ! AND COMMENT THE LOOPS BELOW
1109 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1110 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1112 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1113 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1116 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1118 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1123 ! Substitution for D aminoacids from symmetry.
1126 do j=-nthetyp,nthetyp
1127 do k=-nthetyp,nthetyp
1128 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1130 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1134 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1135 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1136 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1137 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1143 ffthet(llll,lll,ll,i,j,k,iblock)= &
1144 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1145 ffthet(lll,llll,ll,i,j,k,iblock)= &
1146 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1147 ggthet(llll,lll,ll,i,j,k,iblock)= &
1148 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1149 ggthet(lll,llll,ll,i,j,k,iblock)= &
1150 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1159 ! Control printout of the coefficients of virtual-bond-angle potentials
1162 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1167 write (iout,'(//4a)') &
1168 'Type ',onelett(i),onelett(j),onelett(k)
1169 write (iout,'(//a,10x,a)') " l","a[l]"
1170 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1171 write (iout,'(i2,1pe15.5)') &
1172 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1174 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1175 "b",l,"c",l,"d",l,"e",l
1177 write (iout,'(i2,4(1pe15.5))') m,&
1178 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1179 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1183 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1184 "f+",l,"f-",l,"g+",l,"g-",l
1187 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1188 ffthet(n,m,l,i,j,k,iblock),&
1189 ffthet(m,n,l,i,j,k,iblock),&
1190 ggthet(n,m,l,i,j,k,iblock),&
1191 ggthet(m,n,l,i,j,k,iblock)
1201 write (2,*) "Start reading THETA_PDB",ithep_pdb
1203 ! write (2,*) 'i=',i
1204 read (ithep_pdb,*,err=111,end=111) &
1205 a0thet(i),(athet(j,i,1,1),j=1,2),&
1206 (bthet(j,i,1,1),j=1,2)
1207 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1208 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1209 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1210 sigc0(i)=sigc0(i)**2
1213 athet(1,i,1,-1)=athet(1,i,1,1)
1214 athet(2,i,1,-1)=athet(2,i,1,1)
1215 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1216 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1217 athet(1,i,-1,1)=-athet(1,i,1,1)
1218 athet(2,i,-1,1)=-athet(2,i,1,1)
1219 bthet(1,i,-1,1)=bthet(1,i,1,1)
1220 bthet(2,i,-1,1)=bthet(2,i,1,1)
1223 a0thet(i)=a0thet(-i)
1224 athet(1,i,-1,-1)=athet(1,-i,1,1)
1225 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1226 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1227 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1228 athet(1,i,-1,1)=athet(1,-i,1,1)
1229 athet(2,i,-1,1)=-athet(2,-i,1,1)
1230 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1232 athet(1,i,1,-1)=-athet(1,-i,1,1)
1233 athet(2,i,1,-1)=athet(2,-i,1,1)
1234 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1235 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1236 theta0(i)=theta0(-i)
1240 polthet(j,i)=polthet(j,-i)
1243 gthet(j,i)=gthet(j,-i)
1246 write (2,*) "End reading THETA_PDB"
1250 !--------------- Reading theta parameters for nucleic acid-------
1251 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1252 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1253 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1254 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1255 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1256 nthetyp_nucl+1,nthetyp_nucl+1))
1257 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1258 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1259 nthetyp_nucl+1,nthetyp_nucl+1))
1260 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1261 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1262 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1263 nthetyp_nucl+1,nthetyp_nucl+1))
1264 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1265 nthetyp_nucl+1,nthetyp_nucl+1))
1266 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1269 nthetyp_nucl+1,nthetyp_nucl+1))
1270 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1271 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1272 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1273 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1274 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1275 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1277 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1278 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1280 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1282 aa0thet_nucl(:,:,:)=0.0d0
1283 aathet_nucl(:,:,:,:)=0.0d0
1284 bbthet_nucl(:,:,:,:,:)=0.0d0
1285 ccthet_nucl(:,:,:,:,:)=0.0d0
1286 ddthet_nucl(:,:,:,:,:)=0.0d0
1287 eethet_nucl(:,:,:,:,:)=0.0d0
1288 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1289 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1294 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1295 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1296 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1297 read (ithep_nucl,*,end=111,err=111) &
1298 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1299 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1300 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1301 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1302 read (ithep_nucl,*,end=111,err=111) &
1303 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1304 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1305 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1310 !-------------------------------------------
1311 allocate(nlob(ntyp1)) !(ntyp1)
1312 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1313 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1314 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1321 gaussc(:,:,:,:)=0.0D0
1325 ! Read the parameters of the probability distribution/energy expression
1326 ! of the side chains.
1329 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1333 dsc_inv(i)=1.0D0/dsc(i)
1344 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1345 ((blower(k,l,1),l=1,k),k=1,3)
1346 censc(1,1,-i)=censc(1,1,i)
1347 censc(2,1,-i)=censc(2,1,i)
1348 censc(3,1,-i)=-censc(3,1,i)
1350 read (irotam,*,end=112,err=112) bsc(j,i)
1351 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1352 ((blower(k,l,j),l=1,k),k=1,3)
1353 censc(1,j,-i)=censc(1,j,i)
1354 censc(2,j,-i)=censc(2,j,i)
1355 censc(3,j,-i)=-censc(3,j,i)
1356 ! BSC is amplitude of Gaussian
1363 akl=akl+blower(k,m,j)*blower(l,m,j)
1367 if (((k.eq.3).and.(l.ne.3)) &
1368 .or.((l.eq.3).and.(k.ne.3))) then
1369 gaussc(k,l,j,-i)=-akl
1370 gaussc(l,k,j,-i)=-akl
1372 gaussc(k,l,j,-i)=akl
1373 gaussc(l,k,j,-i)=akl
1382 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1385 if (nlobi.gt.0) then
1387 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1388 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1389 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1390 'log h',(bsc(j,i),j=1,nlobi)
1391 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1392 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1394 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1395 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1398 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1399 write (iout,'(a,f10.4,4(16x,f10.4))') &
1400 'Center ',(bsc(j,i),j=1,nlobi)
1401 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1410 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1411 ! added by Urszula Kozlowska 07/11/2007
1413 !el Maximum number of SC local term fitting function coefficiants
1414 !el integer,parameter :: maxsccoef=65
1416 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1419 read (irotam,*,end=112,err=112)
1421 read (irotam,*,end=112,err=112)
1424 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1428 !---------reading nucleic acid parameters for rotamers-------------------
1429 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1430 do i=1,ntyp_molec(2)
1431 read (irotam_nucl,*,end=112,err=112)
1433 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1439 write (iout,*) "Base rotamer parameters"
1440 do i=1,ntyp_molec(2)
1441 write (iout,'(a)') restyp(i,2)
1442 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1447 ! Read the parameters of the probability distribution/energy expression
1448 ! of the side chains.
1450 write (2,*) "Start reading ROTAM_PDB"
1452 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1456 dsc_inv(i)=1.0D0/dsc(i)
1467 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1468 ((blower(k,l,1),l=1,k),k=1,3)
1470 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1471 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1472 ((blower(k,l,j),l=1,k),k=1,3)
1479 akl=akl+blower(k,m,j)*blower(l,m,j)
1489 write (2,*) "End reading ROTAM_PDB"
1495 ! Read torsional parameters in old format
1497 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1499 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1500 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1501 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1503 !el from energy module--------
1504 allocate(v1(nterm_old,ntortyp,ntortyp))
1505 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1506 !el---------------------------
1511 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1517 write (iout,'(/a/)') 'Torsional constants:'
1520 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1521 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1527 ! Read torsional parameters
1529 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1530 read (itorp,*,end=113,err=113) ntortyp
1531 !el from energy module---------
1532 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1533 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1535 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1536 allocate(vlor2(maxlor,ntortyp,ntortyp))
1537 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1538 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1540 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1541 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1542 !el---------------------------
1545 !el---------------------------
1547 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1549 itortyp(i)=-itortyp(-i)
1551 itortyp(ntyp1)=ntortyp
1552 itortyp(-ntyp1)=-ntortyp
1554 write (iout,*) 'ntortyp',ntortyp
1556 do j=-ntortyp+1,ntortyp-1
1557 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1559 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1560 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1563 do k=1,nterm(i,j,iblock)
1564 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1566 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1567 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1568 v0ij=v0ij+si*v1(k,i,j,iblock)
1570 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1571 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1572 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1574 do k=1,nlor(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1576 vlor2(k,i,j),vlor3(k,i,j)
1577 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1580 v0(-i,-j,iblock)=v0ij
1586 write (iout,'(/a/)') 'Torsional constants:'
1588 do i=-ntortyp,ntortyp
1589 do j=-ntortyp,ntortyp
1590 write (iout,*) 'ityp',i,' jtyp',j
1591 write (iout,*) 'Fourier constants'
1592 do k=1,nterm(i,j,iblock)
1593 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1596 write (iout,*) 'Lorenz constants'
1597 do k=1,nlor(i,j,iblock)
1598 write (iout,'(3(1pe15.5))') &
1599 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1605 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1607 ! 6/23/01 Read parameters for double torsionals
1609 !el from energy module------------
1610 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1611 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1612 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1613 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1614 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1615 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1616 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1617 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1618 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1619 !---------------------------------
1623 do j=-ntortyp+1,ntortyp-1
1624 do k=-ntortyp+1,ntortyp-1
1625 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1626 ! write (iout,*) "OK onelett",
1629 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1630 .or. t3.ne.toronelet(k)) then
1631 write (iout,*) "Error in double torsional parameter file",&
1634 call MPI_Finalize(Ierror)
1636 stop "Error in double torsional parameter file"
1638 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1639 ntermd_2(i,j,k,iblock)
1640 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1641 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1642 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1643 ntermd_1(i,j,k,iblock))
1644 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1645 ntermd_1(i,j,k,iblock))
1646 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1647 ntermd_1(i,j,k,iblock))
1648 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1649 ntermd_1(i,j,k,iblock))
1650 ! Martix of D parameters for one dimesional foureir series
1651 do l=1,ntermd_1(i,j,k,iblock)
1652 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1653 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1654 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1655 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1656 ! write(iout,*) "whcodze" ,
1657 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1659 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1660 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1661 v2s(m,l,i,j,k,iblock),&
1662 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1663 ! Martix of D parameters for two dimesional fourier series
1664 do l=1,ntermd_2(i,j,k,iblock)
1666 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1667 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1668 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1669 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1678 write (iout,*) 'Constants for double torsionals'
1681 do j=-ntortyp+1,ntortyp-1
1682 do k=-ntortyp+1,ntortyp-1
1683 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1684 ' nsingle',ntermd_1(i,j,k,iblock),&
1685 ' ndouble',ntermd_2(i,j,k,iblock)
1687 write (iout,*) 'Single angles:'
1688 do l=1,ntermd_1(i,j,k,iblock)
1689 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1690 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1691 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1692 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1695 write (iout,*) 'Pairs of angles:'
1696 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1697 do l=1,ntermd_2(i,j,k,iblock)
1698 write (iout,'(i5,20f10.5)') &
1699 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1702 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1703 do l=1,ntermd_2(i,j,k,iblock)
1704 write (iout,'(i5,20f10.5)') &
1705 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1706 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1715 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1716 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1717 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1718 !el from energy module---------
1719 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1720 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1722 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1723 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1724 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1725 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1727 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1728 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1729 !el---------------------------
1732 !el--------------------
1733 read (itorp_nucl,*,end=113,err=113) &
1734 (itortyp_nucl(i),i=1,ntyp_molec(2))
1735 ! print *,itortyp_nucl(:)
1736 !c write (iout,*) 'ntortyp',ntortyp
1739 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1740 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1743 do k=1,nterm_nucl(i,j)
1744 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1745 v0ij=v0ij+si*v1_nucl(k,i,j)
1748 do k=1,nlor_nucl(i,j)
1749 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1750 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1751 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1757 ! Read of Side-chain backbone correlation parameters
1758 ! Modified 11 May 2012 by Adasko
1761 read (isccor,*,end=119,err=119) nsccortyp
1763 !el from module energy-------------
1764 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1765 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1766 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1767 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1768 !-----------------------------------
1770 !el from module energy-------------
1771 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1773 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1775 isccortyp(i)=-isccortyp(-i)
1777 iscprol=isccortyp(20)
1778 ! write (iout,*) 'ntortyp',ntortyp
1780 !c maxinter is maximum interaction sites
1781 !el from module energy---------
1782 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1783 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1784 -nsccortyp:nsccortyp))
1785 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1786 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1787 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1788 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1789 !-----------------------------------
1793 read (isccor,*,end=119,err=119) &
1794 nterm_sccor(i,j),nlor_sccor(i,j)
1800 nterm_sccor(-i,j)=nterm_sccor(i,j)
1801 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1802 nterm_sccor(i,-j)=nterm_sccor(i,j)
1803 do k=1,nterm_sccor(i,j)
1804 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1806 if (j.eq.iscprol) then
1807 if (i.eq.isccortyp(10)) then
1808 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1809 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1812 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1813 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1814 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1815 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1816 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1817 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1818 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1821 if (i.eq.isccortyp(10)) then
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1823 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1825 if (j.eq.isccortyp(10)) then
1826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1829 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1830 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1831 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1832 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1833 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1834 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1838 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1839 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1840 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1841 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1844 do k=1,nlor_sccor(i,j)
1845 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1846 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1847 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1848 (1+vlor3sccor(k,i,j)**2)
1850 v0sccor(l,i,j)=v0ijsccor
1851 v0sccor(l,-i,j)=v0ijsccor1
1852 v0sccor(l,i,-j)=v0ijsccor2
1853 v0sccor(l,-i,-j)=v0ijsccor3
1859 !el from module energy-------------
1860 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1862 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1863 ! write (iout,*) 'ntortyp',ntortyp
1865 !c maxinter is maximum interaction sites
1866 !el from module energy---------
1867 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1868 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1869 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1870 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1871 !-----------------------------------
1875 read (isccor,*,end=119,err=119) &
1876 nterm_sccor(i,j),nlor_sccor(i,j)
1880 do k=1,nterm_sccor(i,j)
1881 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1883 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1886 do k=1,nlor_sccor(i,j)
1887 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1888 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1889 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1890 (1+vlor3sccor(k,i,j)**2)
1892 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1901 write (iout,'(/a/)') 'Torsional constants:'
1904 write (iout,*) 'ityp',i,' jtyp',j
1905 write (iout,*) 'Fourier constants'
1906 do k=1,nterm_sccor(i,j)
1907 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1909 write (iout,*) 'Lorenz constants'
1910 do k=1,nlor_sccor(i,j)
1911 write (iout,'(3(1pe15.5))') &
1912 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1919 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1920 ! interaction energy of the Gly, Ala, and Pro prototypes.
1925 write (iout,*) "Coefficients of the cumulants"
1927 read (ifourier,*) nloctyp
1928 !write(iout,*) "nloctyp",nloctyp
1929 !el from module energy-------
1930 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1931 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1932 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1933 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1934 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1935 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1936 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1937 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1941 !--------------------------------
1944 read (ifourier,*,end=115,err=115)
1945 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1947 write (iout,*) 'Type',i
1948 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1956 B1tilde(1,i) = buse(3)
1957 B1tilde(2,i) =-buse(5)
1958 B1tilde(1,-i) =-buse(3)
1959 B1tilde(2,-i) =buse(5)
1960 ! buse1tilde(1,i)=0.0d0
1961 ! buse1tilde(2,i)=0.0d0
1981 Ctilde(1,1,i)=buse(7)
1982 Ctilde(1,2,i)=buse(9)
1983 Ctilde(2,1,i)=-buse(9)
1984 Ctilde(2,2,i)=buse(7)
1985 Ctilde(1,1,-i)=buse(7)
1986 Ctilde(1,2,-i)=-buse(9)
1987 Ctilde(2,1,-i)=buse(9)
1988 Ctilde(2,2,-i)=buse(7)
1990 ! Ctilde(1,1,i)=0.0d0
1991 ! Ctilde(1,2,i)=0.0d0
1992 ! Ctilde(2,1,i)=0.0d0
1993 ! Ctilde(2,2,i)=0.0d0
2006 Dtilde(1,1,i)=buse(6)
2007 Dtilde(1,2,i)=buse(8)
2008 Dtilde(2,1,i)=-buse(8)
2009 Dtilde(2,2,i)=buse(6)
2010 Dtilde(1,1,-i)=buse(6)
2011 Dtilde(1,2,-i)=-buse(8)
2012 Dtilde(2,1,-i)=buse(8)
2013 Dtilde(2,2,-i)=buse(6)
2015 ! Dtilde(1,1,i)=0.0d0
2016 ! Dtilde(1,2,i)=0.0d0
2017 ! Dtilde(2,1,i)=0.0d0
2018 ! Dtilde(2,2,i)=0.0d0
2019 EE(1,1,i)= buse(10)+buse(11)
2020 EE(2,2,i)=-buse(10)+buse(11)
2021 EE(2,1,i)= buse(12)-buse(13)
2022 EE(1,2,i)= buse(12)+buse(13)
2023 EE(1,1,-i)= buse(10)+buse(11)
2024 EE(2,2,-i)=-buse(10)+buse(11)
2025 EE(2,1,-i)=-buse(12)+buse(13)
2026 EE(1,2,-i)=-buse(12)-buse(13)
2032 ! ee(2,1,i)=ee(1,2,i)
2036 write (iout,*) 'Type',i
2038 write(iout,*) B1(1,i),B1(2,i)
2040 write(iout,*) B2(1,i),B2(2,i)
2043 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2047 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2051 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2056 ! Read electrostatic-interaction parameters
2061 write (iout,'(/a)') 'Electrostatic interaction constants:'
2062 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2063 'IT','JT','APP','BPP','AEL6','AEL3'
2065 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2066 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2067 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2068 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2073 app (i,j)=epp(i,j)*rri*rri
2074 bpp (i,j)=-2.0D0*epp(i,j)*rri
2075 ael6(i,j)=elpp6(i,j)*4.2D0**6
2076 ael3(i,j)=elpp3(i,j)*4.2D0**3
2078 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2084 ! Read side-chain interaction parameters.
2086 !el from module energy - COMMON.INTERACT-------
2087 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2088 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2089 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2091 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2092 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2093 allocate(epslip(ntyp,ntyp))
2101 !--------------------------------
2103 read (isidep,*,end=117,err=117) ipot,expon
2104 if (ipot.lt.1 .or. ipot.gt.5) then
2105 write (iout,'(2a)') 'Error while reading SC interaction',&
2106 'potential file - unknown potential type.'
2108 call MPI_Finalize(Ierror)
2114 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2115 ', exponents are ',expon,2*expon
2116 ! goto (10,20,30,30,40) ipot
2118 !----------------------- LJ potential ---------------------------------
2120 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2122 (sigma0(i),i=1,ntyp)
2124 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2125 write (iout,'(a/)') 'The epsilon array:'
2126 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2127 write (iout,'(/a)') 'One-body parameters:'
2128 write (iout,'(a,4x,a)') 'residue','sigma'
2129 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2132 !----------------------- LJK potential --------------------------------
2134 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2135 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2136 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2138 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2139 write (iout,'(a/)') 'The epsilon array:'
2140 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2141 write (iout,'(/a)') 'One-body parameters:'
2142 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2143 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2147 !---------------------- GB or BP potential -----------------------------
2151 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2153 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2154 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2155 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2156 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2158 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2161 ! For the GB potential convert sigma'**2 into chi'
2164 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2168 write (iout,'(/a/)') 'Parameters of the BP potential:'
2169 write (iout,'(a/)') 'The epsilon array:'
2170 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2171 write (iout,'(/a)') 'One-body parameters:'
2172 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2174 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2175 chip(i),alp(i),i=1,ntyp)
2178 !--------------------- GBV potential -----------------------------------
2180 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2181 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2182 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2183 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2185 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2186 write (iout,'(a/)') 'The epsilon array:'
2187 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2188 write (iout,'(/a)') 'One-body parameters:'
2189 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2190 's||/s_|_^2',' chip ',' alph '
2191 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2192 sigii(i),chip(i),alp(i),i=1,ntyp)
2195 write(iout,*)"Wrong ipot"
2201 !-----------------------------------------------------------------------
2202 ! Calculate the "working" parameters of SC interactions.
2204 !el from module energy - COMMON.INTERACT-------
2205 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2206 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2207 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2208 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2210 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2223 sc_aa_tube_par(:)=0.0d0
2224 sc_bb_tube_par(:)=0.0d0
2226 !--------------------------------
2231 epslip(i,j)=epslip(j,i)
2236 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2237 sigma(j,i)=sigma(i,j)
2238 rs0(i,j)=dwa16*sigma(i,j)
2242 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2243 'Working parameters of the SC interactions:',&
2244 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2249 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2258 sigeps=dsign(1.0D0,epsij)
2260 aa_aq(i,j)=epsij*rrij*rrij
2261 bb_aq(i,j)=-sigeps*epsij*rrij
2262 aa_aq(j,i)=aa_aq(i,j)
2263 bb_aq(j,i)=bb_aq(i,j)
2264 epsijlip=epslip(i,j)
2265 sigeps=dsign(1.0D0,epsijlip)
2266 epsijlip=dabs(epsijlip)
2267 aa_lip(i,j)=epsijlip*rrij*rrij
2268 bb_lip(i,j)=-sigeps*epsijlip*rrij
2269 aa_lip(j,i)=aa_lip(i,j)
2270 bb_lip(j,i)=bb_lip(i,j)
2271 !C write(iout,*) aa_lip
2273 sigt1sq=sigma0(i)**2
2274 sigt2sq=sigma0(j)**2
2277 ratsig1=sigt2sq/sigt1sq
2278 ratsig2=1.0D0/ratsig1
2279 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2280 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2281 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2285 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2286 sigmaii(i,j)=rsum_max
2287 sigmaii(j,i)=rsum_max
2289 ! sigmaii(i,j)=r0(i,j)
2290 ! sigmaii(j,i)=r0(i,j)
2292 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2293 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2294 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2295 augm(i,j)=epsij*r_augm**(2*expon)
2296 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2303 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2304 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2305 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2310 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2311 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2312 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2313 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2314 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2315 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2316 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2317 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2318 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2320 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2322 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2324 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2325 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2334 read (isidep_nucl,*) ipot_nucl
2335 ! print *,"TU?!",ipot_nucl
2336 if (ipot_nucl.eq.1) then
2337 do i=1,ntyp_molec(2)
2338 do j=i,ntyp_molec(2)
2339 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2340 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2344 do i=1,ntyp_molec(2)
2345 do j=i,ntyp_molec(2)
2346 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2347 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2348 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2352 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2353 do i=1,ntyp_molec(2)
2354 do j=i,ntyp_molec(2)
2355 rrij=sigma_nucl(i,j)
2359 epsij=4*eps_nucl(i,j)
2360 sigeps=dsign(1.0D0,epsij)
2362 aa_nucl(i,j)=epsij*rrij*rrij
2363 bb_nucl(i,j)=-sigeps*epsij*rrij
2364 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2365 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2366 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2367 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2368 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2369 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2372 aa_nucl(i,j)=aa_nucl(j,i)
2373 bb_nucl(i,j)=bb_nucl(j,i)
2374 ael3_nucl(i,j)=ael3_nucl(j,i)
2375 ael6_nucl(i,j)=ael6_nucl(j,i)
2376 ael63_nucl(i,j)=ael63_nucl(j,i)
2377 ael32_nucl(i,j)=ael32_nucl(j,i)
2378 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2379 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2380 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2381 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2382 eps_nucl(i,j)=eps_nucl(j,i)
2383 sigma_nucl(i,j)=sigma_nucl(j,i)
2384 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2388 write(iout,*) "tube param"
2389 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2390 ccavtubpep,dcavtubpep,tubetranenepep
2391 sigmapeptube=sigmapeptube**6
2392 sigeps=dsign(1.0D0,epspeptube)
2393 epspeptube=dabs(epspeptube)
2394 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2395 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2396 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2398 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2399 ccavtub(i),dcavtub(i),tubetranene(i)
2400 sigmasctube=sigmasctube**6
2401 sigeps=dsign(1.0D0,epssctube)
2402 epssctube=dabs(epssctube)
2403 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2404 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2405 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2407 !-----------------READING SC BASE POTENTIALS-----------------------------
2408 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2409 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2410 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2411 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2412 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2413 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2414 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2415 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2416 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2417 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2418 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2419 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2420 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2421 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2422 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2423 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2426 do i=1,ntyp_molec(1)
2427 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2428 write (*,*) "Im in ", i, " ", j
2429 read(isidep_scbase,*) &
2430 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2431 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2432 write(*,*) "eps",eps_scbase(i,j)
2433 read(isidep_scbase,*) &
2434 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2435 chis_scbase(i,j,1),chis_scbase(i,j,2)
2436 read(isidep_scbase,*) &
2437 dhead_scbasei(i,j), &
2438 dhead_scbasej(i,j), &
2439 rborn_scbasei(i,j),rborn_scbasej(i,j)
2440 read(isidep_scbase,*) &
2441 (wdipdip_scbase(k,i,j),k=1,3), &
2442 (wqdip_scbase(k,i,j),k=1,2)
2443 read(isidep_scbase,*) &
2444 alphapol_scbase(i,j), &
2445 epsintab_scbase(i,j)
2448 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2449 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2451 do i=1,ntyp_molec(1)
2452 do j=1,ntyp_molec(2)-1
2453 epsij=eps_scbase(i,j)
2454 rrij=sigma_scbase(i,j)
2459 sigeps=dsign(1.0D0,epsij)
2461 aa_scbase(i,j)=epsij*rrij*rrij
2462 bb_scbase(i,j)=-sigeps*epsij*rrij
2465 !-----------------READING PEP BASE POTENTIALS-------------------
2466 allocate(eps_pepbase(ntyp_molec(2)))
2467 allocate(sigma_pepbase(ntyp_molec(2)))
2468 allocate(chi_pepbase(ntyp_molec(2),2))
2469 allocate(chipp_pepbase(ntyp_molec(2),2))
2470 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2471 allocate(sigmap1_pepbase(ntyp_molec(2)))
2472 allocate(sigmap2_pepbase(ntyp_molec(2)))
2473 allocate(chis_pepbase(ntyp_molec(2),2))
2474 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2477 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2478 write (*,*) "Im in ", i, " ", j
2479 read(isidep_pepbase,*) &
2480 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2481 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2482 write(*,*) "eps",eps_pepbase(j)
2483 read(isidep_pepbase,*) &
2484 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2485 chis_pepbase(j,1),chis_pepbase(j,2)
2486 read(isidep_pepbase,*) &
2487 (wdipdip_pepbase(k,j),k=1,3)
2489 allocate(aa_pepbase(ntyp_molec(2)))
2490 allocate(bb_pepbase(ntyp_molec(2)))
2492 do j=1,ntyp_molec(2)-1
2493 epsij=eps_pepbase(j)
2494 rrij=sigma_pepbase(j)
2499 sigeps=dsign(1.0D0,epsij)
2501 aa_pepbase(j)=epsij*rrij*rrij
2502 bb_pepbase(j)=-sigeps*epsij*rrij
2504 !--------------READING SC PHOSPHATE-------------------------------------
2505 allocate(eps_scpho(ntyp_molec(1)))
2506 allocate(sigma_scpho(ntyp_molec(1)))
2507 allocate(chi_scpho(ntyp_molec(1),2))
2508 allocate(chipp_scpho(ntyp_molec(1),2))
2509 allocate(alphasur_scpho(4,ntyp_molec(1)))
2510 allocate(sigmap1_scpho(ntyp_molec(1)))
2511 allocate(sigmap2_scpho(ntyp_molec(1)))
2512 allocate(chis_scpho(ntyp_molec(1),2))
2513 allocate(wqq_scpho(ntyp_molec(1)))
2514 allocate(wqdip_scpho(2,ntyp_molec(1)))
2515 allocate(alphapol_scpho(ntyp_molec(1)))
2516 allocate(epsintab_scpho(ntyp_molec(1)))
2517 allocate(dhead_scphoi(ntyp_molec(1)))
2518 allocate(rborn_scphoi(ntyp_molec(1)))
2519 allocate(rborn_scphoj(ntyp_molec(1)))
2522 do j=1,ntyp_molec(1) ! without U then we will take T for U
2523 write (*,*) "Im in scpho ", i, " ", j
2524 read(isidep_scpho,*) &
2525 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2526 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2527 write(*,*) "eps",eps_scpho(j)
2528 read(isidep_scpho,*) &
2529 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2530 chis_scpho(j,1),chis_scpho(j,2)
2531 read(isidep_scpho,*) &
2532 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2533 read(isidep_scpho,*) &
2534 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j)
2537 allocate(aa_scpho(ntyp_molec(1)))
2538 allocate(bb_scpho(ntyp_molec(1)))
2540 do j=1,ntyp_molec(1)
2547 sigeps=dsign(1.0D0,epsij)
2549 aa_scpho(j)=epsij*rrij*rrij
2550 bb_scpho(j)=-sigeps*epsij*rrij
2554 read(isidep_peppho,*) &
2555 eps_peppho,sigma_peppho
2556 read(isidep_peppho,*) &
2557 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2558 read(isidep_peppho,*) &
2559 (wqdip_peppho(k),k=1,2)
2567 sigeps=dsign(1.0D0,epsij)
2569 aa_peppho=epsij*rrij*rrij
2570 bb_peppho=-sigeps*epsij*rrij
2573 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2578 ! Define the SC-p interaction constants (hard-coded; old style)
2581 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2583 ! aad(i,1)=0.3D0*4.0D0**12
2584 ! Following line for constants currently implemented
2585 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2586 aad(i,1)=1.5D0*4.0D0**12
2587 ! aad(i,1)=0.17D0*5.6D0**12
2589 ! "Soft" SC-p repulsion
2591 ! Following line for constants currently implemented
2592 ! aad(i,1)=0.3D0*4.0D0**6
2593 ! "Hard" SC-p repulsion
2594 bad(i,1)=3.0D0*4.0D0**6
2595 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2604 ! 8/9/01 Read the SC-p interaction constants from file
2607 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2610 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2611 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2612 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2613 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2617 write (iout,*) "Parameters of SC-p interactions:"
2619 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2620 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2625 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2627 do i=1,ntyp_molec(2)
2628 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2630 do i=1,ntyp_molec(2)
2631 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2632 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2634 r0pp=1.12246204830937298142*5.16158
2640 ! Define the constants of the disulfide bridge
2644 ! Old arbitrary potential - commented out.
2649 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2650 ! energy surface of diethyl disulfide.
2651 ! A. Liwo and U. Kozlowska, 11/24/03
2668 write (iout,'(/a)') "Disulfide bridge parameters:"
2669 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2670 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2671 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2672 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2676 111 write (iout,*) "Error reading bending energy parameters."
2678 112 write (iout,*) "Error reading rotamer energy parameters."
2680 113 write (iout,*) "Error reading torsional energy parameters."
2682 114 write (iout,*) "Error reading double torsional energy parameters."
2684 115 write (iout,*) &
2685 "Error reading cumulant (multibody energy) parameters."
2687 116 write (iout,*) "Error reading electrostatic energy parameters."
2689 117 write (iout,*) "Error reading side chain interaction parameters."
2691 118 write (iout,*) "Error reading SCp interaction parameters."
2693 119 write (iout,*) "Error reading SCCOR parameters"
2696 call MPI_Finalize(Ierror)
2700 end subroutine parmread
2702 !-----------------------------------------------------------------------------
2704 !-----------------------------------------------------------------------------
2705 subroutine printmat(ldim,m,n,iout,key,a)
2708 character(len=3),dimension(n) :: key
2709 real(kind=8),dimension(ldim,n) :: a
2711 integer :: i,j,k,m,iout,nlim
2715 write (iout,1000) (key(k),k=i,nlim)
2717 1000 format (/5x,8(6x,a3))
2718 1020 format (/80(1h-)/)
2720 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2723 1010 format (a3,2x,8(f9.4))
2725 end subroutine printmat
2726 !-----------------------------------------------------------------------------
2728 !-----------------------------------------------------------------------------
2730 ! Read the PDB file and convert the peptide geometry into virtual-chain
2733 use energy_data, only: itype,istype
2737 use control, only: rescode,sugarcode
2738 ! implicit real*8 (a-h,o-z)
2739 ! include 'DIMENSIONS'
2740 ! include 'COMMON.LOCAL'
2741 ! include 'COMMON.VAR'
2742 ! include 'COMMON.CHAIN'
2743 ! include 'COMMON.INTERACT'
2744 ! include 'COMMON.IOUNITS'
2745 ! include 'COMMON.GEO'
2746 ! include 'COMMON.NAMES'
2747 ! include 'COMMON.CONTROL'
2748 ! include 'COMMON.DISTFIT'
2749 ! include 'COMMON.SETUP'
2750 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2752 logical :: lprn=.true.,fail
2753 real(kind=8),dimension(3) :: e1,e2,e3
2754 real(kind=8) :: dcj,efree_temp
2755 character(len=3) :: seq,res
2756 character(len=5) :: atom
2757 character(len=80) :: card
2758 real(kind=8),dimension(3,20) :: sccor
2759 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2760 integer :: isugar,molecprev,firstion
2761 character*1 :: sugar
2763 real(kind=8),dimension(3) :: ccc
2765 integer,dimension(2,maxres/3) :: hfrag_alloc
2766 integer,dimension(4,maxres/3) :: bfrag_alloc
2767 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2768 real(kind=8),dimension(:,:), allocatable :: c_temporary
2769 integer,dimension(:,:) , allocatable :: itype_temporary
2770 integer,dimension(:),allocatable :: istype_temp
2777 ! write (2,*) "UNRES_PDB",unres_pdb
2785 !-----------------------------
2786 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2787 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2790 read (ipdbin,'(a80)',end=10) card
2791 ! write (iout,'(a)') card
2792 if (card(:5).eq.'HELIX') then
2795 read(card(22:25),*) hfrag(1,nhfrag)
2796 read(card(34:37),*) hfrag(2,nhfrag)
2798 if (card(:5).eq.'SHEET') then
2801 read(card(24:26),*) bfrag(1,nbfrag)
2802 read(card(35:37),*) bfrag(2,nbfrag)
2803 !rc----------------------------------------
2804 !rc to be corrected !!!
2805 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2806 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2807 !rc----------------------------------------
2809 if (card(:3).eq.'END') then
2811 else if (card(:3).eq.'TER') then
2816 itype(ires_old,molecule)=ntyp1_molec(molecule)
2817 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2818 nres_molec(molecule)=nres_molec(molecule)+2
2820 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2823 dc(j,ires)=sccor(j,iii)
2826 call sccenter(ires,iii,sccor)
2832 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2833 ! Fish out the ATOM cards.
2834 ! write(iout,*) 'card',card(1:20)
2835 if (index(card(1:4),'ATOM').gt.0) then
2836 read (card(12:16),*) atom
2837 ! write (iout,*) "! ",atom," !",ires
2838 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2839 read (card(23:26),*) ires
2840 read (card(18:20),'(a3)') res
2841 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2842 ! & " ires_old",ires_old
2843 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2844 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2845 if (ires-ishift+ishift1.ne.ires_old) then
2846 ! Calculate the CM of the preceding residue.
2847 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2849 ! write (iout,*) "Calculating sidechain center iii",iii
2852 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2855 call sccenter(ires_old,iii,sccor)
2859 ! Start new residue.
2860 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2863 else if (ibeg.eq.1) then
2864 write (iout,*) "BEG ires",ires
2866 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2869 nres_molec(molecule)=nres_molec(molecule)+1
2871 ires=ires-ishift+ishift1
2873 ! write (iout,*) "ishift",ishift," ires",ires,&
2874 ! " ires_old",ires_old
2876 else if (ibeg.eq.2) then
2878 ishift=-ires_old+ires-1 !!!!!
2879 ishift1=ishift1-1 !!!!!
2880 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2881 ires=ires-ishift+ishift1
2882 ! print *,ires,ishift,ishift1
2886 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2887 ires=ires-ishift+ishift1
2890 ! print *,'atom',ires,atom
2891 if (res.eq.'ACE' .or. res.eq.'NHE') then
2894 if (atom.eq.'CA '.or.atom.eq.'N ') then
2896 itype(ires,molecule)=rescode(ires,res,0,molecule)
2898 ! nres_molec(molecule)=nres_molec(molecule)+1
2901 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2902 ! nres_molec(molecule)=nres_molec(molecule)+1
2903 read (card(19:19),'(a1)') sugar
2904 isugar=sugarcode(sugar,ires)
2905 ! if (ibeg.eq.1) then
2909 ! print *,"ires=",ires,istype(ires)
2915 ires=ires-ishift+ishift1
2917 ! write (iout,*) "ires_old",ires_old," ires",ires
2918 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2921 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2922 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2923 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2924 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2925 ! print *,ires,ishift,ishift1
2926 ! write (iout,*) "backbone ",atom
2928 write (iout,'(2i3,2x,a,3f8.3)') &
2929 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2932 nres_molec(molecule)=nres_molec(molecule)+1
2934 sccor(j,iii)=c(j,ires)
2936 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2937 atom.eq."C2'" .or. atom.eq."C3'" &
2938 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2939 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2940 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2941 ! print *,ires,ishift,ishift1
2945 ! sccor(j,iii)=c(j,ires)
2948 c(j,ires)=c(j,ires)+ccc(j)/5.0
2950 print *,counter,molecule
2951 if (counter.eq.5) then
2953 nres_molec(molecule)=nres_molec(molecule)+1
2956 ! sccor(j,iii)=c(j,ires)
2960 ! print *, "ATOM",atom(1:3)
2961 else if (atom.eq."C5'") then
2962 read (card(19:19),'(a1)') sugar
2963 isugar=sugarcode(sugar,ires)
2968 ! print *,ires,istype(ires)
2971 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2974 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2976 ! write (*,*) card(23:27),ires,itype(ires,1)
2977 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2978 atom.ne.'N' .and. atom.ne.'C' .and. &
2979 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2980 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2981 .and. atom.ne.'P '.and. &
2982 atom(1:1).ne.'H' .and. &
2983 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2985 ! write (iout,*) "sidechain ",atom
2986 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2987 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2988 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2990 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2993 else if ((ions).and.(card(1:6).eq.'HETATM')) then
2994 if (firstion.eq.0) then
2998 dc(j,ires)=sccor(j,iii)
3001 call sccenter(ires,iii,sccor)
3004 read (card(12:16),*) atom
3005 print *,"HETATOM", atom
3006 read (card(18:20),'(a3)') res
3007 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3008 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3009 .or.(atom(1:2).eq.'K ')) &
3012 if (molecule.ne.5) molecprev=molecule
3014 nres_molec(molecule)=nres_molec(molecule)+1
3015 print *,"HERE",nres_molec(molecule)
3016 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
3017 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3021 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3022 if (ires.eq.0) return
3023 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3026 if ((ires_old.ne.ires).and.(molecule.ne.5)) &
3027 nres_molec(molecule)=nres_molec(molecule)-2
3028 ! print *,'I have', nres_molec(:)
3030 do k=1,4 ! ions are without dummy
3031 if (nres_molec(k).eq.0) cycle
3033 ! write (iout,*) i,itype(i,1)
3034 ! if (itype(i,1).eq.ntyp1) then
3035 ! write (iout,*) "dummy",i,itype(i,1)
3037 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3038 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3042 if (itype(i,k).eq.ntyp1_molec(k)) then
3043 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3044 if (itype(i+2,k).eq.0) then
3045 ! print *,"masz sieczke"
3047 if (itype(i+2,j).ne.0) then
3049 itype(i+1,j)=ntyp1_molec(j)
3050 nres_molec(k)=nres_molec(k)-1
3051 nres_molec(j)=nres_molec(j)+1
3057 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3058 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3059 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3060 ! if (unres_pdb) then
3061 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3062 ! print *,i,'tu dochodze'
3063 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3071 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3075 dcj=(c(j,i-2)-c(j,i-3))/2.0
3076 if (dcj.eq.0) dcj=1.23591524223
3081 else !itype(i+1,1).eq.ntyp1
3082 ! if (unres_pdb) then
3083 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3084 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3091 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3095 dcj=(c(j,i+3)-c(j,i+2))/2.0
3096 if (dcj.eq.0) dcj=1.23591524223
3101 endif !itype(i+1,1).eq.ntyp1
3102 endif !itype.eq.ntyp1
3106 ! Calculate the CM of the last side chain.
3110 dc(j,ires)=sccor(j,iii)
3113 call sccenter(ires,iii,sccor)
3119 ! print *,"molecule",molecule
3120 if ((itype(nres,1).ne.10)) then
3122 if (molecule.eq.5) molecule=molecprev
3123 itype(nres,molecule)=ntyp1_molec(molecule)
3124 nres_molec(molecule)=nres_molec(molecule)+1
3125 ! if (unres_pdb) then
3126 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3127 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3134 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3138 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3139 c(j,nres)=c(j,nres-1)+dcj
3140 c(j,2*nres)=c(j,nres)
3144 ! print *,'I have',nres, nres_molec(:)
3146 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3148 if (nres.ne.nres0) then
3149 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3151 stop "Error nres value in WHAM input"
3154 !---------------------------------
3155 !el reallocate tables
3158 ! hfrag_alloc(j,i)=hfrag(j,i)
3161 ! bfrag_alloc(j,i)=bfrag(j,i)
3167 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3168 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3169 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3170 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3174 ! hfrag(j,i)=hfrag_alloc(j,i)
3179 ! bfrag(j,i)=bfrag_alloc(j,i)
3182 !el end reallocate tables
3183 !---------------------------------
3191 c(j,2*nres)=c(j,nres)
3194 if (itype(1,1).eq.ntyp1) then
3198 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3199 call refsys(2,3,4,e1,e2,e3,fail)
3206 c(j,1)=c(j,2)-1.9d0*e2(j)
3210 dcj=(c(j,4)-c(j,3))/2.0
3216 ! First lets assign correct dummy to correct type of chain
3218 if (itype(1,1).eq.ntyp1) then
3219 if (itype(2,1).eq.0) then
3221 if (itype(2,j).ne.0) then
3223 itype(1,j)=ntyp1_molec(j)
3224 nres_molec(1)=nres_molec(1)-1
3225 nres_molec(j)=nres_molec(j)+1
3232 print *,'I have',nres, nres_molec(:)
3234 ! Copy the coordinates to reference coordinates
3240 ! Calculate internal coordinates.
3242 write (iout,'(/a)') &
3243 "Cartesian coordinates of the reference structure"
3244 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3245 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3247 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3248 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3249 (c(j,ires+nres),j=1,3)
3252 ! znamy już nres wiec mozna alokowac tablice
3253 ! Calculate internal coordinates.
3254 if(me.eq.king.or..not.out1file)then
3255 write (iout,'(a)') &
3256 "Backbone and SC coordinates as read from the PDB"
3258 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3259 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3260 (c(j,nres+ires),j=1,3)
3263 ! NOW LETS ROCK! SORTING
3264 allocate(c_temporary(3,2*nres))
3265 allocate(itype_temporary(nres,5))
3266 allocate(molnum(nres+1))
3267 allocate(istype_temp(nres))
3268 itype_temporary(:,:)=0
3272 if (itype(i,k).ne.0) then
3274 c_temporary(j,seqalingbegin)=c(j,i)
3275 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3278 itype_temporary(seqalingbegin,k)=itype(i,k)
3279 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3280 istype_temp(seqalingbegin)=istype(i)
3281 molnum(seqalingbegin)=k
3282 seqalingbegin=seqalingbegin+1
3288 c(j,i)=c_temporary(j,i)
3293 itype(i,k)=itype_temporary(i,k)
3294 istype(i)=istype_temp(i)
3297 if (itype(1,1).eq.ntyp1) then
3301 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3302 call refsys(2,3,4,e1,e2,e3,fail)
3309 c(j,1)=c(j,2)-1.9d0*e2(j)
3313 dcj=(c(j,4)-c(j,3))/2.0
3321 write (iout,'(/a)') &
3322 "Cartesian coordinates of the reference structure after sorting"
3323 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3324 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3326 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3327 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3328 (c(j,ires+nres),j=1,3)
3332 ! print *,seqalingbegin,nres
3333 if(.not.allocated(vbld)) then
3334 allocate(vbld(2*nres))
3339 if(.not.allocated(vbld_inv)) then
3340 allocate(vbld_inv(2*nres))
3346 if(.not.allocated(theta)) then
3347 allocate(theta(nres+2))
3351 if(.not.allocated(phi)) allocate(phi(nres+2))
3352 if(.not.allocated(alph)) allocate(alph(nres+2))
3353 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3354 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3355 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3356 if(.not.allocated(costtab)) allocate(costtab(nres))
3357 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3358 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3359 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3360 if(.not.allocated(xxref)) allocate(xxref(nres))
3361 if(.not.allocated(yyref)) allocate(yyref(nres))
3362 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3363 if(.not.allocated(dc_norm)) then
3364 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3365 allocate(dc_norm(3,0:2*nres+2))
3369 call int_from_cart(.true.,.false.)
3370 call sc_loc_geom(.false.)
3372 thetaref(i)=theta(i)
3382 dc(j,i)=c(j,i+1)-c(j,i)
3383 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3388 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3389 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3391 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3395 ! Copy the coordinates to reference coordinates
3396 ! Splits to single chain if occurs
3400 ! cref(j,i,cou)=c(j,i)
3404 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3405 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3406 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3407 !-----------------------------
3411 write (iout,*) "symetr", symetr
3414 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3416 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3419 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3424 cref(j,i,cou)=c(j,i)
3425 cref(j,i+nres,cou)=c(j,i+nres)
3427 chain_rep(j,lll,kkk)=c(j,i)
3428 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3432 write (iout,*) chain_length
3433 if (chain_length.eq.0) chain_length=nres
3435 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3436 chain_rep(j,chain_length+nres,symetr) &
3437 =chain_rep(j,chain_length+nres,1)
3440 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3442 ! do kkk=1,chain_length
3443 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3447 ! makes copy of chains
3448 write (iout,*) "symetr", symetr
3453 if (symetr.gt.1) then
3460 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3466 ! write (iout,*) i,icha
3467 do lll=1,chain_length
3469 if (cou.le.nres) then
3471 kupa=mod(lll,chain_length)
3472 iprzes=(kkk-1)*chain_length+lll
3473 if (kupa.eq.0) kupa=chain_length
3474 ! write (iout,*) "kupa", kupa
3475 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3476 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3483 !-koniec robienia kopii
3486 write (iout,*) "nowa struktura", nperm
3488 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3490 cref(3,i,kkk),cref(1,nres+i,kkk),&
3491 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3493 100 format (//' alpha-carbon coordinates ',&
3494 ' centroid coordinates'/ &
3495 ' ', 6X,'X',11X,'Y',11X,'Z', &
3496 10X,'X',11X,'Y',11X,'Z')
3497 110 format (a,'(',i3,')',6f12.5)
3503 bfrag(i,j)=bfrag(i,j)-ishift
3509 hfrag(i,j)=hfrag(i,j)-ishift
3515 end subroutine readpdb
3516 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3517 !-----------------------------------------------------------------------------
3519 !-----------------------------------------------------------------------------
3520 subroutine read_control
3534 use random, only: random_init
3535 ! implicit real*8 (a-h,o-z)
3536 ! include 'DIMENSIONS'
3538 use prng, only:prng_restart
3540 logical :: OKRandom!, prng_restart
3543 ! include 'COMMON.IOUNITS'
3544 ! include 'COMMON.TIME1'
3545 ! include 'COMMON.THREAD'
3546 ! include 'COMMON.SBRIDGE'
3547 ! include 'COMMON.CONTROL'
3548 ! include 'COMMON.MCM'
3549 ! include 'COMMON.MAP'
3550 ! include 'COMMON.HEADER'
3551 ! include 'COMMON.CSA'
3552 ! include 'COMMON.CHAIN'
3553 ! include 'COMMON.MUCA'
3554 ! include 'COMMON.MD'
3555 ! include 'COMMON.FFIELD'
3556 ! include 'COMMON.INTERACT'
3557 ! include 'COMMON.SETUP'
3558 !el integer :: KDIAG,ICORFL,IXDR
3559 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3560 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3561 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3562 ! character(len=80) :: ucase
3563 character(len=640) :: controlcard
3565 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3571 read (INP,'(a)') titel
3572 call card_concat(controlcard,.true.)
3573 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3574 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3575 call reada(controlcard,'SEED',seed,0.0D0)
3576 call random_init(seed)
3577 ! Set up the time limit (caution! The time must be input in minutes!)
3578 read_cart=index(controlcard,'READ_CART').gt.0
3579 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3580 call readi(controlcard,'SYM',symetr,1)
3581 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3582 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3583 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3584 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3585 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3586 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3587 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3588 call reada(controlcard,'DRMS',drms,0.1D0)
3589 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3590 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3591 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3592 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3593 write (iout,'(a,f10.1)')'DRMS = ',drms
3594 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3595 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3597 call readi(controlcard,'NZ_START',nz_start,0)
3598 call readi(controlcard,'NZ_END',nz_end,0)
3599 call readi(controlcard,'IZ_SC',iz_sc,0)
3600 timlim=60.0D0*timlim
3601 safety = 60.0d0*safety
3604 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3605 !C SHIELD keyword sets if the shielding effect of side-chains is used
3606 !C 0 denots no shielding is used all peptide are equally despite the
3607 !C solvent accesible area
3608 !C 1 the newly introduced function
3609 !C 2 reseved for further possible developement
3610 call readi(controlcard,'SHIELD',shield_mode,0)
3611 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3612 write(iout,*) "shield_mode",shield_mode
3613 !C Varibles set size of box
3614 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3615 protein=index(controlcard,"PROTEIN").gt.0
3616 ions=index(controlcard,"IONS").gt.0
3617 nucleic=index(controlcard,"NUCLEIC").gt.0
3618 write (iout,*) "with_theta_constr ",with_theta_constr
3619 AFMlog=(index(controlcard,'AFM'))
3620 selfguide=(index(controlcard,'SELFGUIDE'))
3621 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3622 call readi(controlcard,'GENCONSTR',genconstr,0)
3623 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3624 call reada(controlcard,'BOXY',boxysize,100.0d0)
3625 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3626 call readi(controlcard,'TUBEMOD',tubemode,0)
3627 write (iout,*) TUBEmode,"TUBEMODE"
3628 if (TUBEmode.gt.0) then
3629 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3630 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3631 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3632 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3633 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3634 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3635 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3636 buftubebot=bordtubebot+tubebufthick
3637 buftubetop=bordtubetop-tubebufthick
3640 ! CUTOFFF ON ELECTROSTATICS
3641 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3642 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3643 write(iout,*) "R_CUT_ELE=",r_cut_ele
3644 ! Lipidic parameters
3645 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3646 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3647 if (lipthick.gt.0.0d0) then
3648 bordliptop=(boxzsize+lipthick)/2.0
3649 bordlipbot=bordliptop-lipthick
3650 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3651 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3652 buflipbot=bordlipbot+lipbufthick
3653 bufliptop=bordliptop-lipbufthick
3654 if ((lipbufthick*2.0d0).gt.lipthick) &
3655 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3656 endif !lipthick.gt.0
3657 write(iout,*) "bordliptop=",bordliptop
3658 write(iout,*) "bordlipbot=",bordlipbot
3659 write(iout,*) "bufliptop=",bufliptop
3660 write(iout,*) "buflipbot=",buflipbot
3661 write (iout,*) "SHIELD MODE",shield_mode
3663 !C-------------------------
3664 minim=(index(controlcard,'MINIMIZE').gt.0)
3665 dccart=(index(controlcard,'CART').gt.0)
3666 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3667 overlapsc=.not.overlapsc
3668 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3669 searchsc=.not.searchsc
3670 sideadd=(index(controlcard,'SIDEADD').gt.0)
3671 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3672 outpdb=(index(controlcard,'PDBOUT').gt.0)
3673 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3674 pdbref=(index(controlcard,'PDBREF').gt.0)
3675 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3676 indpdb=index(controlcard,'PDBSTART')
3677 extconf=(index(controlcard,'EXTCONF').gt.0)
3678 call readi(controlcard,'IPRINT',iprint,0)
3679 call readi(controlcard,'MAXGEN',maxgen,10000)
3680 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3681 call readi(controlcard,"KDIAG",kdiag,0)
3682 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3683 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3684 write (iout,*) "RESCALE_MODE",rescale_mode
3685 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3686 if (index(controlcard,'REGULAR').gt.0.0D0) then
3687 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3691 if (index(controlcard,'CHECKGRAD').gt.0) then
3693 if (index(controlcard,'CART').gt.0) then
3695 elseif (index(controlcard,'CARINT').gt.0) then
3700 elseif (index(controlcard,'THREAD').gt.0) then
3702 call readi(controlcard,'THREAD',nthread,0)
3703 if (nthread.gt.0) then
3704 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3707 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3708 stop 'Error termination in Read_Control.'
3710 else if (index(controlcard,'MCMA').gt.0) then
3712 else if (index(controlcard,'MCEE').gt.0) then
3714 else if (index(controlcard,'MULTCONF').gt.0) then
3716 else if (index(controlcard,'MAP').gt.0) then
3718 call readi(controlcard,'MAP',nmap,0)
3719 else if (index(controlcard,'CSA').gt.0) then
3721 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3723 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3726 !fcm else if (index(controlcard,'MCMF').gt.0) then
3728 else if (index(controlcard,'SOFTREG').gt.0) then
3730 else if (index(controlcard,'CHECK_BOND').gt.0) then
3732 else if (index(controlcard,'TEST').gt.0) then
3734 else if (index(controlcard,'MD').gt.0) then
3736 else if (index(controlcard,'RE ').gt.0) then
3740 lmuca=index(controlcard,'MUCA').gt.0
3741 call readi(controlcard,'MUCADYN',mucadyn,0)
3742 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3743 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3745 write (iout,*) 'MUCADYN=',mucadyn
3746 write (iout,*) 'MUCASMOOTH=',muca_smooth
3749 iscode=index(controlcard,'ONE_LETTER')
3750 indphi=index(controlcard,'PHI')
3751 indback=index(controlcard,'BACK')
3752 iranconf=index(controlcard,'RAND_CONF')
3753 i2ndstr=index(controlcard,'USE_SEC_PRED')
3754 gradout=index(controlcard,'GRADOUT').gt.0
3755 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3756 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3757 if (me.eq.king .or. .not.out1file ) &
3758 write (iout,*) "DISTCHAINMAX",distchainmax
3760 if(me.eq.king.or..not.out1file) &
3761 write (iout,'(2a)') diagmeth(kdiag),&
3762 ' routine used to diagonalize matrices.'
3763 if (shield_mode.gt.0) then
3765 !C VSolvSphere the volume of solving sphere
3767 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3768 !C there will be no distinction between proline peptide group and normal peptide
3769 !C group in case of shielding parameters
3770 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3771 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3772 write (iout,*) VSolvSphere,VSolvSphere_div
3773 !C long axis of side chain
3775 long_r_sidechain(i)=vbldsc0(1,i)
3776 short_r_sidechain(i)=sigma0(i)
3777 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3783 end subroutine read_control
3784 !-----------------------------------------------------------------------------
3785 subroutine read_REMDpar
3787 ! Read REMD settings
3794 use control_data, only:out1file
3795 ! implicit real*8 (a-h,o-z)
3796 ! include 'DIMENSIONS'
3797 ! include 'COMMON.IOUNITS'
3798 ! include 'COMMON.TIME1'
3799 ! include 'COMMON.MD'
3802 !el include 'COMMON.LANGEVIN'
3804 !el include 'COMMON.LANGEVIN.lang0'
3806 ! include 'COMMON.INTERACT'
3807 ! include 'COMMON.NAMES'
3808 ! include 'COMMON.GEO'
3809 ! include 'COMMON.REMD'
3810 ! include 'COMMON.CONTROL'
3811 ! include 'COMMON.SETUP'
3812 ! character(len=80) :: ucase
3813 character(len=320) :: controlcard
3814 character(len=3200) :: controlcard1
3815 integer :: iremd_m_total
3818 ! real(kind=8) :: var,ene
3820 if(me.eq.king.or..not.out1file) &
3821 write (iout,*) "REMD setup"
3823 call card_concat(controlcard,.true.)
3824 call readi(controlcard,"NREP",nrep,3)
3825 call readi(controlcard,"NSTEX",nstex,1000)
3826 call reada(controlcard,"RETMIN",retmin,10.0d0)
3827 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3828 mremdsync=(index(controlcard,'SYNC').gt.0)
3829 call readi(controlcard,"NSYN",i_sync_step,100)
3830 restart1file=(index(controlcard,'REST1FILE').gt.0)
3831 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3832 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3833 if(max_cache_traj_use.gt.max_cache_traj) &
3834 max_cache_traj_use=max_cache_traj
3835 if(me.eq.king.or..not.out1file) then
3836 !d if (traj1file) then
3837 !rc caching is in testing - NTWX is not ignored
3838 !d write (iout,*) "NTWX value is ignored"
3839 !d write (iout,*) " trajectory is stored to one file by master"
3840 !d write (iout,*) " before exchange at NSTEX intervals"
3842 write (iout,*) "NREP= ",nrep
3843 write (iout,*) "NSTEX= ",nstex
3844 write (iout,*) "SYNC= ",mremdsync
3845 write (iout,*) "NSYN= ",i_sync_step
3846 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3849 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3850 if (index(controlcard,'TLIST').gt.0) then
3852 call card_concat(controlcard1,.true.)
3853 read(controlcard1,*) (remd_t(i),i=1,nrep)
3854 if(me.eq.king.or..not.out1file) &
3855 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3858 if (index(controlcard,'MLIST').gt.0) then
3860 call card_concat(controlcard1,.true.)
3861 read(controlcard1,*) (remd_m(i),i=1,nrep)
3862 if(me.eq.king.or..not.out1file) then
3863 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3866 iremd_m_total=iremd_m_total+remd_m(i)
3868 write (iout,*) 'Total number of replicas ',iremd_m_total
3871 if(me.eq.king.or..not.out1file) &
3872 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3874 end subroutine read_REMDpar
3875 !-----------------------------------------------------------------------------
3876 subroutine read_MDpar
3880 use control_data, only: r_cut,rlamb,out1file
3882 use geometry_data, only: pi
3884 ! implicit real*8 (a-h,o-z)
3885 ! include 'DIMENSIONS'
3886 ! include 'COMMON.IOUNITS'
3887 ! include 'COMMON.TIME1'
3888 ! include 'COMMON.MD'
3891 !el include 'COMMON.LANGEVIN'
3893 !el include 'COMMON.LANGEVIN.lang0'
3895 ! include 'COMMON.INTERACT'
3896 ! include 'COMMON.NAMES'
3897 ! include 'COMMON.GEO'
3898 ! include 'COMMON.SETUP'
3899 ! include 'COMMON.CONTROL'
3900 ! include 'COMMON.SPLITELE'
3901 ! character(len=80) :: ucase
3902 character(len=320) :: controlcard
3907 call card_concat(controlcard,.true.)
3908 call readi(controlcard,"NSTEP",n_timestep,1000000)
3909 call readi(controlcard,"NTWE",ntwe,100)
3910 call readi(controlcard,"NTWX",ntwx,1000)
3911 call reada(controlcard,"DT",d_time,1.0d-1)
3912 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3913 call reada(controlcard,"DAMAX",damax,1.0d1)
3914 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3915 call readi(controlcard,"LANG",lang,0)
3916 RESPA = index(controlcard,"RESPA") .gt. 0
3917 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3918 ntime_split0=ntime_split
3919 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3920 ntime_split0=ntime_split
3921 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3922 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3923 rest = index(controlcard,"REST").gt.0
3924 tbf = index(controlcard,"TBF").gt.0
3925 usampl = index(controlcard,"USAMPL").gt.0
3926 mdpdb = index(controlcard,"MDPDB").gt.0
3927 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3928 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3929 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3930 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3931 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3932 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3933 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3934 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3935 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3936 large = index(controlcard,"LARGE").gt.0
3937 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3938 rattle = index(controlcard,"RATTLE").gt.0
3939 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3945 if(me.eq.king.or..not.out1file) then
3947 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3949 write (iout,'(a)') "The units are:"
3950 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3951 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3952 " acceleration: angstrom/(48.9 fs)**2"
3953 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3955 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3956 write (iout,'(a60,f10.5,a)') &
3957 "Initial time step of numerical integration:",d_time,&
3959 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3961 write (iout,'(2a,i4,a)') &
3962 "A-MTS algorithm used; initial time step for fast-varying",&
3963 " short-range forces split into",ntime_split," steps."
3964 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3965 r_cut," lambda",rlamb
3967 write (iout,'(2a,f10.5)') &
3968 "Maximum acceleration threshold to reduce the time step",&
3969 "/increase split number:",damax
3970 write (iout,'(2a,f10.5)') &
3971 "Maximum predicted energy drift to reduce the timestep",&
3972 "/increase split number:",edriftmax
3973 write (iout,'(a60,f10.5)') &
3974 "Maximum velocity threshold to reduce velocities:",dvmax
3975 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3976 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3977 if (rattle) write (iout,'(a60)') &
3978 "Rattle algorithm used to constrain the virtual bonds"
3982 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3983 call reada(controlcard,"RWAT",rwat,1.4d0)
3984 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3985 surfarea=index(controlcard,"SURFAREA").gt.0
3986 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3987 if(me.eq.king.or..not.out1file)then
3988 write (iout,'(/a,$)') "Langevin dynamics calculation"
3990 write (iout,'(a/)') &
3991 " with direct integration of Langevin equations"
3992 else if (lang.eq.2) then
3993 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3994 else if (lang.eq.3) then
3995 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3996 else if (lang.eq.4) then
3997 write (iout,'(a/)') " in overdamped mode"
3999 write (iout,'(//a,i5)') &
4000 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4003 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4004 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4005 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4006 write (iout,'(a60,f10.5)') &
4007 "Scaling factor of the friction forces:",scal_fric
4008 if (surfarea) write (iout,'(2a,i10,a)') &
4009 "Friction coefficients will be scaled by solvent-accessible",&
4010 " surface area every",reset_fricmat," steps."
4012 ! Calculate friction coefficients and bounds of stochastic forces
4013 eta=6*pi*cPoise*etawat
4014 if(me.eq.king.or..not.out1file) &
4015 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4018 do j=1,5 !types of molecules
4019 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4020 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4022 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4023 do j=1,5 !types of molecules
4025 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4026 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4030 if(me.eq.king.or..not.out1file)then
4031 write (iout,'(/2a/)') &
4032 "Radii of site types and friction coefficients and std's of",&
4033 " stochastic forces of fully exposed sites"
4034 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4036 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4037 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4041 if(me.eq.king.or..not.out1file)then
4042 write (iout,'(a)') "Berendsen bath calculation"
4043 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4044 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4046 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4047 count_reset_moment," steps"
4049 write (iout,'(a,i10,a)') &
4050 "Velocities will be reset at random every",count_reset_vel,&
4054 if(me.eq.king.or..not.out1file) &
4055 write (iout,'(a31)') "Microcanonical mode calculation"
4057 if(me.eq.king.or..not.out1file)then
4058 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4060 write(iout,*) "MD running with constraints."
4061 write(iout,*) "Equilibration time ", eq_time, " mtus."
4062 write(iout,*) "Constraining ", nfrag," fragments."
4063 write(iout,*) "Length of each fragment, weight and q0:"
4065 write (iout,*) "Set of restraints #",iset
4067 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4068 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4070 write(iout,*) "constraints between ", npair, "fragments."
4071 write(iout,*) "constraint pairs, weights and q0:"
4073 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4074 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4076 write(iout,*) "angle constraints within ", nfrag_back,&
4077 "backbone fragments."
4078 write(iout,*) "fragment, weights:"
4080 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4081 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4082 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4085 iset=mod(kolor,nset)+1
4088 if(me.eq.king.or..not.out1file) &
4089 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4091 end subroutine read_MDpar
4092 !-----------------------------------------------------------------------------
4096 ! implicit real*8 (a-h,o-z)
4097 ! include 'DIMENSIONS'
4098 ! include 'COMMON.MAP'
4099 ! include 'COMMON.IOUNITS'
4100 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4101 character(len=80) :: mapcard !,ucase
4104 ! real(kind=8) :: var,ene
4107 read (inp,'(a)') mapcard
4108 mapcard=ucase(mapcard)
4109 if (index(mapcard,'PHI').gt.0) then
4111 else if (index(mapcard,'THE').gt.0) then
4113 else if (index(mapcard,'ALP').gt.0) then
4115 else if (index(mapcard,'OME').gt.0) then
4118 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4119 stop 'Error - illegal variable spec in MAP card.'
4121 call readi (mapcard,'RES1',res1(imap),0)
4122 call readi (mapcard,'RES2',res2(imap),0)
4123 if (res1(imap).eq.0) then
4124 res1(imap)=res2(imap)
4125 else if (res2(imap).eq.0) then
4126 res2(imap)=res1(imap)
4128 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4129 write (iout,'(a)') &
4130 'Error - illegal definition of variable group in MAP.'
4131 stop 'Error - illegal definition of variable group in MAP.'
4133 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4134 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4135 call readi(mapcard,'NSTEP',nstep(imap),0)
4136 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4137 write (iout,'(a)') &
4138 'Illegal boundary and/or step size specification in MAP.'
4139 stop 'Illegal boundary and/or step size specification in MAP.'
4143 end subroutine map_read
4144 !-----------------------------------------------------------------------------
4147 use control_data, only: vdisulf
4149 ! implicit real*8 (a-h,o-z)
4150 ! include 'DIMENSIONS'
4151 ! include 'COMMON.IOUNITS'
4152 ! include 'COMMON.GEO'
4153 ! include 'COMMON.CSA'
4154 ! include 'COMMON.BANK'
4155 ! include 'COMMON.CONTROL'
4156 ! character(len=80) :: ucase
4157 character(len=620) :: mcmcard
4159 ! integer :: ntf,ik,iw_pdb
4160 ! real(kind=8) :: var,ene
4162 call card_concat(mcmcard,.true.)
4164 call readi(mcmcard,'NCONF',nconf,50)
4165 call readi(mcmcard,'NADD',nadd,0)
4166 call readi(mcmcard,'JSTART',jstart,1)
4167 call readi(mcmcard,'JEND',jend,1)
4168 call readi(mcmcard,'NSTMAX',nstmax,500000)
4169 call readi(mcmcard,'N0',n0,1)
4170 call readi(mcmcard,'N1',n1,6)
4171 call readi(mcmcard,'N2',n2,4)
4172 call readi(mcmcard,'N3',n3,0)
4173 call readi(mcmcard,'N4',n4,0)
4174 call readi(mcmcard,'N5',n5,0)
4175 call readi(mcmcard,'N6',n6,10)
4176 call readi(mcmcard,'N7',n7,0)
4177 call readi(mcmcard,'N8',n8,0)
4178 call readi(mcmcard,'N9',n9,0)
4179 call readi(mcmcard,'N14',n14,0)
4180 call readi(mcmcard,'N15',n15,0)
4181 call readi(mcmcard,'N16',n16,0)
4182 call readi(mcmcard,'N17',n17,0)
4183 call readi(mcmcard,'N18',n18,0)
4185 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4187 call readi(mcmcard,'NDIFF',ndiff,2)
4188 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4189 call readi(mcmcard,'IS1',is1,1)
4190 call readi(mcmcard,'IS2',is2,8)
4191 call readi(mcmcard,'NRAN0',nran0,4)
4192 call readi(mcmcard,'NRAN1',nran1,2)
4193 call readi(mcmcard,'IRR',irr,1)
4194 call readi(mcmcard,'NSEED',nseed,20)
4195 call readi(mcmcard,'NTOTAL',ntotal,10000)
4196 call reada(mcmcard,'CUT1',cut1,2.0d0)
4197 call reada(mcmcard,'CUT2',cut2,5.0d0)
4198 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4199 call readi(mcmcard,'ICMAX',icmax,3)
4200 call readi(mcmcard,'IRESTART',irestart,0)
4201 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4204 call reada(mcmcard,'DELE',dele,20.0d0)
4205 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4206 call readi(mcmcard,'IREF',iref,0)
4207 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4208 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4209 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4210 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4211 write (iout,*) "NCONF_IN",nconf_in
4213 end subroutine csaread
4214 !-----------------------------------------------------------------------------
4218 use control_data, only: MaxMoveType
4221 ! implicit real*8 (a-h,o-z)
4222 ! include 'DIMENSIONS'
4223 ! include 'COMMON.MCM'
4224 ! include 'COMMON.MCE'
4225 ! include 'COMMON.IOUNITS'
4226 ! character(len=80) :: ucase
4227 character(len=320) :: mcmcard
4230 ! real(kind=8) :: var,ene
4232 call card_concat(mcmcard,.true.)
4233 call readi(mcmcard,'MAXACC',maxacc,100)
4234 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4235 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4236 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4237 call readi(mcmcard,'MAXREPM',maxrepm,200)
4238 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4239 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4240 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4241 call reada(mcmcard,'E_UP',e_up,5.0D0)
4242 call reada(mcmcard,'DELTE',delte,0.1D0)
4243 call readi(mcmcard,'NSWEEP',nsweep,5)
4244 call readi(mcmcard,'NSTEPH',nsteph,0)
4245 call readi(mcmcard,'NSTEPC',nstepc,0)
4246 call reada(mcmcard,'TMIN',tmin,298.0D0)
4247 call reada(mcmcard,'TMAX',tmax,298.0D0)
4248 call readi(mcmcard,'NWINDOW',nwindow,0)
4249 call readi(mcmcard,'PRINT_MC',print_mc,0)
4250 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4251 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4252 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4253 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4254 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4255 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4256 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4257 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4258 if (nwindow.gt.0) then
4259 allocate(winstart(nwindow)) !!el (maxres)
4260 allocate(winend(nwindow)) !!el
4261 allocate(winlen(nwindow)) !!el
4262 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4264 winlen(i)=winend(i)-winstart(i)+1
4267 if (tmax.lt.tmin) tmax=tmin
4268 if (tmax.eq.tmin) then
4272 if (nstepc.gt.0 .and. nsteph.gt.0) then
4273 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4274 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4276 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4277 ! Probabilities of different move types
4278 sumpro_type(0)=0.0D0
4279 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4280 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4281 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4282 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4283 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4284 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4285 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4287 print *,'i',i,' sumprotype',sumpro_type(i)
4288 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4289 print *,'i',i,' sumprotype',sumpro_type(i)
4292 end subroutine mcmread
4293 !-----------------------------------------------------------------------------
4294 subroutine read_minim
4297 ! implicit real*8 (a-h,o-z)
4298 ! include 'DIMENSIONS'
4299 ! include 'COMMON.MINIM'
4300 ! include 'COMMON.IOUNITS'
4301 ! character(len=80) :: ucase
4302 character(len=320) :: minimcard
4304 ! integer :: ntf,ik,iw_pdb
4305 ! real(kind=8) :: var,ene
4307 call card_concat(minimcard,.true.)
4308 call readi(minimcard,'MAXMIN',maxmin,2000)
4309 call readi(minimcard,'MAXFUN',maxfun,5000)
4310 call readi(minimcard,'MINMIN',minmin,maxmin)
4311 call readi(minimcard,'MINFUN',minfun,maxmin)
4312 call reada(minimcard,'TOLF',tolf,1.0D-2)
4313 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4314 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4315 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4316 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4317 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4318 'Options in energy minimization:'
4319 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4320 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4321 'MinMin:',MinMin,' MinFun:',MinFun,&
4322 ' TolF:',TolF,' RTolF:',RTolF
4324 end subroutine read_minim
4325 !-----------------------------------------------------------------------------
4326 subroutine openunits
4328 use MD_data, only: usampl
4331 use control_data, only:out1file
4332 use control, only: getenv_loc
4333 ! implicit real*8 (a-h,o-z)
4334 ! include 'DIMENSIONS'
4337 character(len=16) :: form,nodename
4338 integer :: nodelen,ierror,npos
4340 ! include 'COMMON.SETUP'
4341 ! include 'COMMON.IOUNITS'
4342 ! include 'COMMON.MD'
4343 ! include 'COMMON.CONTROL'
4344 integer :: lenpre,lenpot,lentmp !,ilen
4346 character(len=3) :: out1file_text !,ucase
4347 character(len=3) :: ll
4350 ! integer :: ntf,ik,iw_pdb
4351 ! real(kind=8) :: var,ene
4353 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4354 call getenv_loc("PREFIX",prefix)
4356 call getenv_loc("POT",pot)
4357 call getenv_loc("DIRTMP",tmpdir)
4358 call getenv_loc("CURDIR",curdir)
4359 call getenv_loc("OUT1FILE",out1file_text)
4360 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4361 out1file_text=ucase(out1file_text)
4362 if (out1file_text(1:1).eq."Y") then
4365 out1file=fg_rank.gt.0
4370 if (lentmp.gt.0) then
4371 write (*,'(80(1h!))')
4372 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4373 write (*,'(80(1h!))')
4374 write (*,*)"All output files will be on node /tmp directory."
4376 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4377 if (me.eq.king) then
4378 write (*,*) "The master node is ",nodename
4379 else if (fg_rank.eq.0) then
4380 write (*,*) "I am the CG slave node ",nodename
4382 write (*,*) "I am the FG slave node ",nodename
4385 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4386 lenpre = lentmp+lenpre+1
4388 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4389 ! Get the names and open the input files
4390 #if defined(WINIFL) || defined(WINPGI)
4391 open(1,file=pref_orig(:ilen(pref_orig))// &
4392 '.inp',status='old',readonly,shared)
4393 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4394 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4395 ! Get parameter filenames and open the parameter files.
4396 call getenv_loc('BONDPAR',bondname)
4397 open (ibond,file=bondname,status='old',readonly,shared)
4398 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4399 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4400 call getenv_loc('THETPAR',thetname)
4401 open (ithep,file=thetname,status='old',readonly,shared)
4402 call getenv_loc('ROTPAR',rotname)
4403 open (irotam,file=rotname,status='old',readonly,shared)
4404 call getenv_loc('TORPAR',torname)
4405 open (itorp,file=torname,status='old',readonly,shared)
4406 call getenv_loc('TORDPAR',tordname)
4407 open (itordp,file=tordname,status='old',readonly,shared)
4408 call getenv_loc('FOURIER',fouriername)
4409 open (ifourier,file=fouriername,status='old',readonly,shared)
4410 call getenv_loc('ELEPAR',elename)
4411 open (ielep,file=elename,status='old',readonly,shared)
4412 call getenv_loc('SIDEPAR',sidename)
4413 open (isidep,file=sidename,status='old',readonly,shared)
4415 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4416 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4417 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4418 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4419 call getenv_loc('TORPAR_NUCL',torname_nucl)
4420 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4421 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4422 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4423 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4424 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4427 #elif (defined CRAY) || (defined AIX)
4428 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4430 ! print *,"Processor",myrank," opened file 1"
4431 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4432 ! print *,"Processor",myrank," opened file 9"
4433 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4434 ! Get parameter filenames and open the parameter files.
4435 call getenv_loc('BONDPAR',bondname)
4436 open (ibond,file=bondname,status='old',action='read')
4437 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4438 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4440 ! print *,"Processor",myrank," opened file IBOND"
4441 call getenv_loc('THETPAR',thetname)
4442 open (ithep,file=thetname,status='old',action='read')
4443 ! print *,"Processor",myrank," opened file ITHEP"
4444 call getenv_loc('ROTPAR',rotname)
4445 open (irotam,file=rotname,status='old',action='read')
4446 ! print *,"Processor",myrank," opened file IROTAM"
4447 call getenv_loc('TORPAR',torname)
4448 open (itorp,file=torname,status='old',action='read')
4449 ! print *,"Processor",myrank," opened file ITORP"
4450 call getenv_loc('TORDPAR',tordname)
4451 open (itordp,file=tordname,status='old',action='read')
4452 ! print *,"Processor",myrank," opened file ITORDP"
4453 call getenv_loc('SCCORPAR',sccorname)
4454 open (isccor,file=sccorname,status='old',action='read')
4455 ! print *,"Processor",myrank," opened file ISCCOR"
4456 call getenv_loc('FOURIER',fouriername)
4457 open (ifourier,file=fouriername,status='old',action='read')
4458 ! print *,"Processor",myrank," opened file IFOURIER"
4459 call getenv_loc('ELEPAR',elename)
4460 open (ielep,file=elename,status='old',action='read')
4461 ! print *,"Processor",myrank," opened file IELEP"
4462 call getenv_loc('SIDEPAR',sidename)
4463 open (isidep,file=sidename,status='old',action='read')
4465 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4466 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4467 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4468 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4469 call getenv_loc('TORPAR_NUCL',torname_nucl)
4470 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4471 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4472 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4473 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4474 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4476 call getenv_loc('LIPTRANPAR',liptranname)
4477 open (iliptranpar,file=liptranname,status='old',action='read')
4478 call getenv_loc('TUBEPAR',tubename)
4479 open (itube,file=tubename,status='old',action='read')
4480 call getenv_loc('IONPAR',ionname)
4481 open (iion,file=ionname,status='old',action='read')
4483 ! print *,"Processor",myrank," opened file ISIDEP"
4484 ! print *,"Processor",myrank," opened parameter files"
4486 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4487 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4488 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4489 ! Get parameter filenames and open the parameter files.
4490 call getenv_loc('BONDPAR',bondname)
4491 open (ibond,file=bondname,status='old')
4492 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4493 open (ibond_nucl,file=bondname_nucl,status='old')
4495 call getenv_loc('THETPAR',thetname)
4496 open (ithep,file=thetname,status='old')
4497 call getenv_loc('ROTPAR',rotname)
4498 open (irotam,file=rotname,status='old')
4499 call getenv_loc('TORPAR',torname)
4500 open (itorp,file=torname,status='old')
4501 call getenv_loc('TORDPAR',tordname)
4502 open (itordp,file=tordname,status='old')
4503 call getenv_loc('SCCORPAR',sccorname)
4504 open (isccor,file=sccorname,status='old')
4505 call getenv_loc('FOURIER',fouriername)
4506 open (ifourier,file=fouriername,status='old')
4507 call getenv_loc('ELEPAR',elename)
4508 open (ielep,file=elename,status='old')
4509 call getenv_loc('SIDEPAR',sidename)
4510 open (isidep,file=sidename,status='old')
4512 open (ithep_nucl,file=thetname_nucl,status='old')
4513 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4514 open (irotam_nucl,file=rotname_nucl,status='old')
4515 call getenv_loc('TORPAR_NUCL',torname_nucl)
4516 open (itorp_nucl,file=torname_nucl,status='old')
4517 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4518 open (itordp_nucl,file=tordname_nucl,status='old')
4519 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4520 open (isidep_nucl,file=sidename_nucl,status='old')
4522 call getenv_loc('LIPTRANPAR',liptranname)
4523 open (iliptranpar,file=liptranname,status='old')
4524 call getenv_loc('TUBEPAR',tubename)
4525 open (itube,file=tubename,status='old')
4526 call getenv_loc('IONPAR',ionname)
4527 open (iion,file=ionname,status='old')
4529 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4531 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4532 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4533 ! Get parameter filenames and open the parameter files.
4534 call getenv_loc('BONDPAR',bondname)
4535 open (ibond,file=bondname,status='old',action='read')
4536 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4537 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4538 call getenv_loc('THETPAR',thetname)
4539 open (ithep,file=thetname,status='old',action='read')
4540 call getenv_loc('ROTPAR',rotname)
4541 open (irotam,file=rotname,status='old',action='read')
4542 call getenv_loc('TORPAR',torname)
4543 open (itorp,file=torname,status='old',action='read')
4544 call getenv_loc('TORDPAR',tordname)
4545 open (itordp,file=tordname,status='old',action='read')
4546 call getenv_loc('SCCORPAR',sccorname)
4547 open (isccor,file=sccorname,status='old',action='read')
4549 call getenv_loc('THETPARPDB',thetname_pdb)
4550 print *,"thetname_pdb ",thetname_pdb
4551 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4552 print *,ithep_pdb," opened"
4554 call getenv_loc('FOURIER',fouriername)
4555 open (ifourier,file=fouriername,status='old',readonly)
4556 call getenv_loc('ELEPAR',elename)
4557 open (ielep,file=elename,status='old',readonly)
4558 call getenv_loc('SIDEPAR',sidename)
4559 open (isidep,file=sidename,status='old',readonly)
4561 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4562 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4563 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4564 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4565 call getenv_loc('TORPAR_NUCL',torname_nucl)
4566 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4567 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4568 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4569 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4570 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4571 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4572 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4573 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4574 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4575 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4576 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4577 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4578 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4581 call getenv_loc('LIPTRANPAR',liptranname)
4582 open (iliptranpar,file=liptranname,status='old',action='read')
4583 call getenv_loc('TUBEPAR',tubename)
4584 open (itube,file=tubename,status='old',action='read')
4585 call getenv_loc('IONPAR',ionname)
4586 open (iion,file=ionname,status='old',action='read')
4589 call getenv_loc('ROTPARPDB',rotname_pdb)
4590 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4593 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4594 #if defined(WINIFL) || defined(WINPGI)
4595 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4596 #elif (defined CRAY) || (defined AIX)
4597 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4599 open (iscpp_nucl,file=scpname_nucl,status='old')
4601 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4606 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4607 ! Use -DOLDSCP to use hard-coded constants instead.
4609 call getenv_loc('SCPPAR',scpname)
4610 #if defined(WINIFL) || defined(WINPGI)
4611 open (iscpp,file=scpname,status='old',readonly,shared)
4612 #elif (defined CRAY) || (defined AIX)
4613 open (iscpp,file=scpname,status='old',action='read')
4615 open (iscpp,file=scpname,status='old')
4617 open (iscpp,file=scpname,status='old',action='read')
4620 call getenv_loc('PATTERN',patname)
4621 #if defined(WINIFL) || defined(WINPGI)
4622 open (icbase,file=patname,status='old',readonly,shared)
4623 #elif (defined CRAY) || (defined AIX)
4624 open (icbase,file=patname,status='old',action='read')
4626 open (icbase,file=patname,status='old')
4628 open (icbase,file=patname,status='old',action='read')
4631 ! Open output file only for CG processes
4632 ! print *,"Processor",myrank," fg_rank",fg_rank
4633 if (fg_rank.eq.0) then
4635 if (nodes.eq.1) then
4638 npos = dlog10(dfloat(nodes-1))+1
4640 if (npos.lt.3) npos=3
4641 write (liczba,'(i1)') npos
4642 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4644 write (liczba,form) me
4645 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4646 liczba(:ilen(liczba))
4647 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4649 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4651 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4652 liczba(:ilen(liczba))//'.mol2'
4653 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4654 liczba(:ilen(liczba))//'.stat'
4656 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4657 //liczba(:ilen(liczba))//'.stat')
4658 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4661 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4662 liczba(:ilen(liczba))//'.const'
4667 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4668 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4669 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4670 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4671 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4673 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4675 rest2name=prefix(:ilen(prefix))//'.rst'
4677 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4680 #if defined(AIX) || defined(PGI)
4681 if (me.eq.king .or. .not. out1file) &
4682 open(iout,file=outname,status='unknown')
4684 if (fg_rank.gt.0) then
4685 write (liczba,'(i3.3)') myrank/nfgtasks
4686 write (ll,'(bz,i3.3)') fg_rank
4687 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4692 open(igeom,file=intname,status='unknown',position='append')
4693 open(ipdb,file=pdbname,status='unknown')
4694 open(imol2,file=mol2name,status='unknown')
4695 open(istat,file=statname,status='unknown',position='append')
4697 !1out open(iout,file=outname,status='unknown')
4700 if (me.eq.king .or. .not.out1file) &
4701 open(iout,file=outname,status='unknown')
4703 if (fg_rank.gt.0) then
4704 write (liczba,'(i3.3)') myrank/nfgtasks
4705 write (ll,'(bz,i3.3)') fg_rank
4706 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4711 open(igeom,file=intname,status='unknown',access='append')
4712 open(ipdb,file=pdbname,status='unknown')
4713 open(imol2,file=mol2name,status='unknown')
4714 open(istat,file=statname,status='unknown',access='append')
4716 !1out open(iout,file=outname,status='unknown')
4719 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4720 csa_seed=prefix(:lenpre)//'.CSA.seed'
4721 csa_history=prefix(:lenpre)//'.CSA.history'
4722 csa_bank=prefix(:lenpre)//'.CSA.bank'
4723 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4724 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4725 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4726 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4727 csa_int=prefix(:lenpre)//'.int'
4728 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4729 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4730 csa_in=prefix(:lenpre)//'.CSA.in'
4731 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4734 write (iout,'(80(1h-))')
4735 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4736 write (iout,'(80(1h-))')
4737 write (iout,*) "Input file : ",&
4738 pref_orig(:ilen(pref_orig))//'.inp'
4739 write (iout,*) "Output file : ",&
4740 outname(:ilen(outname))
4742 write (iout,*) "Sidechain potential file : ",&
4743 sidename(:ilen(sidename))
4745 write (iout,*) "SCp potential file : ",&
4746 scpname(:ilen(scpname))
4748 write (iout,*) "Electrostatic potential file : ",&
4749 elename(:ilen(elename))
4750 write (iout,*) "Cumulant coefficient file : ",&
4751 fouriername(:ilen(fouriername))
4752 write (iout,*) "Torsional parameter file : ",&
4753 torname(:ilen(torname))
4754 write (iout,*) "Double torsional parameter file : ",&
4755 tordname(:ilen(tordname))
4756 write (iout,*) "SCCOR parameter file : ",&
4757 sccorname(:ilen(sccorname))
4758 write (iout,*) "Bond & inertia constant file : ",&
4759 bondname(:ilen(bondname))
4760 write (iout,*) "Bending parameter file : ",&
4761 thetname(:ilen(thetname))
4762 write (iout,*) "Rotamer parameter file : ",&
4763 rotname(:ilen(rotname))
4766 write (iout,*) "Thetpdb parameter file : ",&
4767 thetname_pdb(:ilen(thetname_pdb))
4770 write (iout,*) "Threading database : ",&
4771 patname(:ilen(patname))
4773 write (iout,*)" DIRTMP : ",&
4775 write (iout,'(80(1h-))')
4778 end subroutine openunits
4779 !-----------------------------------------------------------------------------
4782 use geometry_data, only: nres,dc
4784 ! implicit real*8 (a-h,o-z)
4785 ! include 'DIMENSIONS'
4786 ! include 'COMMON.CHAIN'
4787 ! include 'COMMON.IOUNITS'
4788 ! include 'COMMON.MD'
4791 ! real(kind=8) :: var,ene
4793 open(irest2,file=rest2name,status='unknown')
4794 read(irest2,*) totT,EK,potE,totE,t_bath
4797 ! AL 4/17/17: Now reading d_t(0,:) too
4799 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4802 ! AL 4/17/17: Now reading d_c(0,:) too
4804 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4807 read (irest2,*) iset
4811 end subroutine readrst
4812 !-----------------------------------------------------------------------------
4813 subroutine read_fragments
4817 use control_data, only:out1file
4820 ! implicit real*8 (a-h,o-z)
4821 ! include 'DIMENSIONS'
4825 ! include 'COMMON.SETUP'
4826 ! include 'COMMON.CHAIN'
4827 ! include 'COMMON.IOUNITS'
4828 ! include 'COMMON.MD'
4829 ! include 'COMMON.CONTROL'
4832 ! real(kind=8) :: var,ene
4834 read(inp,*) nset,nfrag,npair,nfrag_back
4836 !el from module energy
4837 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4838 if(.not.allocated(wfrag_back)) then
4839 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4840 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4842 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4843 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4845 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4846 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4849 if(me.eq.king.or..not.out1file) &
4850 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4851 " nfrag_back",nfrag_back
4853 read(inp,*) mset(iset)
4855 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4857 if(me.eq.king.or..not.out1file) &
4858 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4859 ifrag(2,i,iset), qinfrag(i,iset)
4862 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4864 if(me.eq.king.or..not.out1file) &
4865 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4866 ipair(2,i,iset), qinpair(i,iset)
4869 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4870 wfrag_back(3,i,iset),&
4871 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4872 if(me.eq.king.or..not.out1file) &
4873 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4874 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4878 end subroutine read_fragments
4879 !-----------------------------------------------------------------------------
4881 !-----------------------------------------------------------------------------
4885 ! implicit real*8 (a-h,o-z)
4886 ! include 'DIMENSIONS'
4887 ! include 'COMMON.CSA'
4888 ! include 'COMMON.BANK'
4889 ! include 'COMMON.IOUNITS'
4891 ! integer :: ntf,ik,iw_pdb
4892 ! real(kind=8) :: var,ene
4894 open(icsa_in,file=csa_in,status="old",err=100)
4895 read(icsa_in,*) nconf
4896 read(icsa_in,*) jstart,jend
4897 read(icsa_in,*) nstmax
4898 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4899 read(icsa_in,*) nran0,nran1,irr
4900 read(icsa_in,*) nseed
4901 read(icsa_in,*) ntotal,cut1,cut2
4902 read(icsa_in,*) estop
4903 read(icsa_in,*) icmax,irestart
4904 read(icsa_in,*) ntbankm,dele,difcut
4905 read(icsa_in,*) iref,rmscut,pnccut
4906 read(icsa_in,*) ndiff
4913 end subroutine csa_read
4914 !-----------------------------------------------------------------------------
4915 subroutine initial_write
4918 ! implicit real*8 (a-h,o-z)
4919 ! include 'DIMENSIONS'
4920 ! include 'COMMON.CSA'
4921 ! include 'COMMON.BANK'
4922 ! include 'COMMON.IOUNITS'
4924 ! integer :: ntf,ik,iw_pdb
4925 ! real(kind=8) :: var,ene
4927 open(icsa_seed,file=csa_seed,status="unknown")
4928 write(icsa_seed,*) "seed"
4930 #if defined(AIX) || defined(PGI)
4931 open(icsa_history,file=csa_history,status="unknown",&
4934 open(icsa_history,file=csa_history,status="unknown",&
4937 write(icsa_history,*) nconf
4938 write(icsa_history,*) jstart,jend
4939 write(icsa_history,*) nstmax
4940 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4941 write(icsa_history,*) nran0,nran1,irr
4942 write(icsa_history,*) nseed
4943 write(icsa_history,*) ntotal,cut1,cut2
4944 write(icsa_history,*) estop
4945 write(icsa_history,*) icmax,irestart
4946 write(icsa_history,*) ntbankm,dele,difcut
4947 write(icsa_history,*) iref,rmscut,pnccut
4948 write(icsa_history,*) ndiff
4950 write(icsa_history,*)
4953 open(icsa_bank1,file=csa_bank1,status="unknown")
4954 write(icsa_bank1,*) 0
4958 end subroutine initial_write
4959 !-----------------------------------------------------------------------------
4960 subroutine restart_write
4963 ! implicit real*8 (a-h,o-z)
4964 ! include 'DIMENSIONS'
4965 ! include 'COMMON.IOUNITS'
4966 ! include 'COMMON.CSA'
4967 ! include 'COMMON.BANK'
4969 ! integer :: ntf,ik,iw_pdb
4970 ! real(kind=8) :: var,ene
4972 #if defined(AIX) || defined(PGI)
4973 open(icsa_history,file=csa_history,position="append")
4975 open(icsa_history,file=csa_history,access="append")
4977 write(icsa_history,*)
4978 write(icsa_history,*) "This is restart"
4979 write(icsa_history,*)
4980 write(icsa_history,*) nconf
4981 write(icsa_history,*) jstart,jend
4982 write(icsa_history,*) nstmax
4983 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4984 write(icsa_history,*) nran0,nran1,irr
4985 write(icsa_history,*) nseed
4986 write(icsa_history,*) ntotal,cut1,cut2
4987 write(icsa_history,*) estop
4988 write(icsa_history,*) icmax,irestart
4989 write(icsa_history,*) ntbankm,dele,difcut
4990 write(icsa_history,*) iref,rmscut,pnccut
4991 write(icsa_history,*) ndiff
4992 write(icsa_history,*)
4993 write(icsa_history,*) "irestart is: ", irestart
4995 write(icsa_history,*)
4999 end subroutine restart_write
5000 !-----------------------------------------------------------------------------
5002 !-----------------------------------------------------------------------------
5003 subroutine write_pdb(npdb,titelloc,ee)
5005 ! implicit real*8 (a-h,o-z)
5006 ! include 'DIMENSIONS'
5007 ! include 'COMMON.IOUNITS'
5008 character(len=50) :: titelloc1
5009 character*(*) :: titelloc
5010 character(len=3) :: zahl
5011 character(len=5) :: liczba5
5013 integer :: npdb !,ilen
5017 ! real(kind=8) :: var,ene
5021 if (npdb.lt.1000) then
5022 call numstr(npdb,zahl)
5023 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5025 if (npdb.lt.10000) then
5026 write(liczba5,'(i1,i4)') 0,npdb
5028 write(liczba5,'(i5)') npdb
5030 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5032 call pdbout(ee,titelloc1,ipdb)
5035 end subroutine write_pdb
5036 !-----------------------------------------------------------------------------
5038 !-----------------------------------------------------------------------------
5039 subroutine write_thread_summary
5040 ! Thread the sequence through a database of known structures
5041 use control_data, only: refstr
5043 use energy_data, only: n_ene_comp
5045 ! implicit real*8 (a-h,o-z)
5046 ! include 'DIMENSIONS'
5048 use MPI_data !include 'COMMON.INFO'
5051 ! include 'COMMON.CONTROL'
5052 ! include 'COMMON.CHAIN'
5053 ! include 'COMMON.DBASE'
5054 ! include 'COMMON.INTERACT'
5055 ! include 'COMMON.VAR'
5056 ! include 'COMMON.THREAD'
5057 ! include 'COMMON.FFIELD'
5058 ! include 'COMMON.SBRIDGE'
5059 ! include 'COMMON.HEADER'
5060 ! include 'COMMON.NAMES'
5061 ! include 'COMMON.IOUNITS'
5062 ! include 'COMMON.TIME1'
5064 integer,dimension(maxthread) :: ip
5065 real(kind=8),dimension(0:n_ene) :: energia
5067 integer :: i,j,ii,jj,ipj,ik,kk,ist
5068 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5070 write (iout,'(30x,a/)') &
5071 ' *********** Summary threading statistics ************'
5072 write (iout,'(a)') 'Initial energies:'
5073 write (iout,'(a4,2x,a12,14a14,3a8)') &
5074 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5075 'RMSnat','NatCONT','NNCONT','RMS'
5076 ! Energy sort patterns
5081 enet=ener(n_ene-1,ip(i))
5084 if (ener(n_ene-1,ip(j)).lt.enet) then
5086 enet=ener(n_ene-1,ip(j))
5098 ist=nres_base(2,ii)+ipatt(2,i)
5100 energia(i)=ener0(kk,i)
5102 etot=ener0(n_ene_comp+1,i)
5103 rmsnat=ener0(n_ene_comp+2,i)
5104 rms=ener0(n_ene_comp+3,i)
5105 frac=ener0(n_ene_comp+4,i)
5106 frac_nn=ener0(n_ene_comp+5,i)
5109 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5110 i,str_nam(ii),ist+1,&
5111 (energia(print_order(kk)),kk=1,nprint_ene),&
5112 etot,rmsnat,frac,frac_nn,rms
5114 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5115 i,str_nam(ii),ist+1,&
5116 (energia(print_order(kk)),kk=1,nprint_ene),etot
5119 write (iout,'(//a)') 'Final energies:'
5120 write (iout,'(a4,2x,a12,17a14,3a8)') &
5121 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5122 'RMSnat','NatCONT','NNCONT','RMS'
5126 ist=nres_base(2,ii)+ipatt(2,i)
5128 energia(kk)=ener(kk,ik)
5130 etot=ener(n_ene_comp+1,i)
5131 rmsnat=ener(n_ene_comp+2,i)
5132 rms=ener(n_ene_comp+3,i)
5133 frac=ener(n_ene_comp+4,i)
5134 frac_nn=ener(n_ene_comp+5,i)
5135 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5136 i,str_nam(ii),ist+1,&
5137 (energia(print_order(kk)),kk=1,nprint_ene),&
5138 etot,rmsnat,frac,frac_nn,rms
5140 write (iout,'(/a/)') 'IEXAM array:'
5141 write (iout,'(i5)') nexcl
5143 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5145 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5146 'Max. time for threading step ',max_time_for_thread,&
5147 'Average time for threading step: ',ave_time_for_thread
5149 end subroutine write_thread_summary
5150 !-----------------------------------------------------------------------------
5151 subroutine write_stat_thread(ithread,ipattern,ist)
5153 use energy_data, only: n_ene_comp
5155 ! implicit real*8 (a-h,o-z)
5156 ! include "DIMENSIONS"
5157 ! include "COMMON.CONTROL"
5158 ! include "COMMON.IOUNITS"
5159 ! include "COMMON.THREAD"
5160 ! include "COMMON.FFIELD"
5161 ! include "COMMON.DBASE"
5162 ! include "COMMON.NAMES"
5163 real(kind=8),dimension(0:n_ene) :: energia
5165 integer :: ithread,ipattern,ist,i
5166 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5168 #if defined(AIX) || defined(PGI)
5169 open(istat,file=statname,position='append')
5171 open(istat,file=statname,access='append')
5174 energia(i)=ener(i,ithread)
5176 etot=ener(n_ene_comp+1,ithread)
5177 rmsnat=ener(n_ene_comp+2,ithread)
5178 rms=ener(n_ene_comp+3,ithread)
5179 frac=ener(n_ene_comp+4,ithread)
5180 frac_nn=ener(n_ene_comp+5,ithread)
5181 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5182 ithread,str_nam(ipattern),ist+1,&
5183 (energia(print_order(i)),i=1,nprint_ene),&
5184 etot,rmsnat,frac,frac_nn,rms
5187 end subroutine write_stat_thread
5188 !-----------------------------------------------------------------------------
5190 !-----------------------------------------------------------------------------
5191 end module io_config