8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 sigmasctube,krad2,ract
749 integer :: ichir1,ichir2,ijunk,irdiff
751 character*80 temp1,mychar
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
817 allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
818 allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
834 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
837 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
838 ! dsc(i) = vbldsc0_nucl(1,i)
842 ! dsc_inv(i)=1.0D0/dsc(i)
845 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
846 ! do i=1,ntyp_molec(2)
847 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
848 ! aksc_nucl(j,i),abond0_nucl(j,i),&
849 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
850 ! dsc(i) = vbldsc0(1,i)
854 ! dsc_inv(i)=1.0D0/dsc(i)
859 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
860 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
862 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
864 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
865 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
867 write (iout,'(13x,3f10.5)') &
868 vbldsc0(j,i),aksc(j,i),abond0(j,i)
875 if (.not.allocated(ichargecat)) &
876 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
878 if (oldion.eq.1) then
880 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
881 print *,msc(i,5),restok(i,5)
885 read (iion,*) ncatprotparm
886 allocate(catprm(ncatprotparm,4))
888 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
892 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
896 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
899 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
900 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
904 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
905 restyp(i,5),"-",restyp(j,2)
908 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
909 !----------------------------------------------------
910 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
911 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
912 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
913 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
914 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
915 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
925 allocate(liptranene(ntyp))
926 !C reading lipid parameters
927 write (iout,*) "iliptranpar",iliptranpar
929 read(iliptranpar,*) pepliptran
932 read(iliptranpar,*) liptranene(i)
933 print *,liptranene(i)
936 ! water parmaters entalphy
937 allocate(awaterenta(0:400))
938 allocate(bwaterenta(0:400))
939 allocate(cwaterenta(0:400))
940 allocate(dwaterenta(0:400))
941 allocate(awaterentro(0:400))
942 allocate(bwaterentro(0:400))
943 allocate(cwaterentro(0:400))
944 allocate(dwaterentro(0:400))
946 read(iwaterwater,*) mychar
947 read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
948 cwaterenta(0),dwaterenta(0)
950 read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
951 cwaterenta(i),dwaterenta(i)
952 irdiff=int((ract-2.06d0)*50.0d0)+1
953 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
955 ! water parmaters entrophy
956 read(iwaterwater,*) mychar
957 read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
958 cwaterentro(0),dwaterentro(0)
960 read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
961 cwaterentro(i),dwaterentro(i)
962 irdiff=int((ract-2.06d0)*50.0d0)+1
963 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
969 ! Read the parameters of the probability distribution/energy expression
970 ! of the virtual-bond valence angles theta
973 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
974 (bthet(j,i,1,1),j=1,2)
975 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
976 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
977 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
981 athet(1,i,1,-1)=athet(1,i,1,1)
982 athet(2,i,1,-1)=athet(2,i,1,1)
983 bthet(1,i,1,-1)=-bthet(1,i,1,1)
984 bthet(2,i,1,-1)=-bthet(2,i,1,1)
985 athet(1,i,-1,1)=-athet(1,i,1,1)
986 athet(2,i,-1,1)=-athet(2,i,1,1)
987 bthet(1,i,-1,1)=bthet(1,i,1,1)
988 bthet(2,i,-1,1)=bthet(2,i,1,1)
992 athet(1,i,-1,-1)=athet(1,-i,1,1)
993 athet(2,i,-1,-1)=-athet(2,-i,1,1)
994 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
995 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
996 athet(1,i,-1,1)=athet(1,-i,1,1)
997 athet(2,i,-1,1)=-athet(2,-i,1,1)
998 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
999 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1000 athet(1,i,1,-1)=-athet(1,-i,1,1)
1001 athet(2,i,1,-1)=athet(2,-i,1,1)
1002 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1003 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1004 theta0(i)=theta0(-i)
1008 polthet(j,i)=polthet(j,-i)
1011 gthet(j,i)=gthet(j,-i)
1017 if (.not.LaTeX) then
1018 write (iout,'(a)') &
1019 'Parameters of the virtual-bond valence angles:'
1020 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
1021 ' ATHETA0 ',' A1 ',' A2 ',&
1024 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1025 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
1027 write (iout,'(/a/9x,5a/79(1h-))') &
1028 'Parameters of the expression for sigma(theta_c):',&
1029 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
1030 ' ALPH3 ',' SIGMA0C '
1032 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1033 (polthet(j,i),j=0,3),sigc0(i)
1035 write (iout,'(/a/9x,5a/79(1h-))') &
1036 'Parameters of the second gaussian:',&
1037 ' THETA0 ',' SIGMA0 ',' G1 ',&
1040 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1041 sig0(i),(gthet(j,i),j=1,3)
1044 write (iout,'(a)') &
1045 'Parameters of the virtual-bond valence angles:'
1046 write (iout,'(/a/9x,5a/79(1h-))') &
1047 'Coefficients of expansion',&
1048 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1049 ' b1*10^1 ',' b2*10^1 '
1051 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1052 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1053 (10*bthet(j,i,1,1),j=1,2)
1055 write (iout,'(/a/9x,5a/79(1h-))') &
1056 'Parameters of the expression for sigma(theta_c):',&
1057 ' alpha0 ',' alph1 ',' alph2 ',&
1058 ' alhp3 ',' sigma0c '
1060 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1061 (polthet(j,i),j=0,3),sigc0(i)
1063 write (iout,'(/a/9x,5a/79(1h-))') &
1064 'Parameters of the second gaussian:',&
1065 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1068 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1069 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1075 ! Read the parameters of Utheta determined from ab initio surfaces
1076 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1078 IF (tor_mode.eq.0) THEN
1079 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1080 ntheterm3,nsingle,ndouble
1081 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1083 !----------------------------------------------------
1084 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1085 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1086 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1087 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1088 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1089 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1090 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1091 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1092 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1093 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1094 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1095 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1096 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1097 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1098 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1099 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1100 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1101 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1102 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1103 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1104 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1105 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1106 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1107 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1109 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1111 ithetyp(i)=-ithetyp(-i)
1114 aa0thet(:,:,:,:)=0.0d0
1115 aathet(:,:,:,:,:)=0.0d0
1116 bbthet(:,:,:,:,:,:)=0.0d0
1117 ccthet(:,:,:,:,:,:)=0.0d0
1118 ddthet(:,:,:,:,:,:)=0.0d0
1119 eethet(:,:,:,:,:,:)=0.0d0
1120 ffthet(:,:,:,:,:,:,:)=0.0d0
1121 ggthet(:,:,:,:,:,:,:)=0.0d0
1123 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1125 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1126 ! VAR:1=non-glicyne non-proline 2=proline
1127 ! VAR:negative values for D-aminoacid
1129 do j=-nthetyp,nthetyp
1130 do k=-nthetyp,nthetyp
1131 read (ithep,'(6a)',end=111,err=111) res1
1132 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1133 ! VAR: aa0thet is variable describing the average value of Foureir
1134 ! VAR: expansion series
1135 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1136 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1137 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1138 read (ithep,*,end=111,err=111) &
1139 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1140 read (ithep,*,end=111,err=111) &
1141 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1142 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1143 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1144 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1146 read (ithep,*,end=111,err=111) &
1147 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1148 ffthet(lll,llll,ll,i,j,k,iblock),&
1149 ggthet(llll,lll,ll,i,j,k,iblock),&
1150 ggthet(lll,llll,ll,i,j,k,iblock),&
1151 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1156 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1157 ! coefficients of theta-and-gamma-dependent terms are zero.
1158 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1159 ! RECOMENTDED AFTER VERSION 3.3)
1163 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1164 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1166 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1167 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1170 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1172 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1175 ! AND COMMENT THE LOOPS BELOW
1179 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1180 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1182 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1183 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1186 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1188 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1193 ! Substitution for D aminoacids from symmetry.
1196 do j=-nthetyp,nthetyp
1197 do k=-nthetyp,nthetyp
1198 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1200 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1204 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1205 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1206 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1207 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1213 ffthet(llll,lll,ll,i,j,k,iblock)= &
1214 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1215 ffthet(lll,llll,ll,i,j,k,iblock)= &
1216 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1217 ggthet(llll,lll,ll,i,j,k,iblock)= &
1218 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1219 ggthet(lll,llll,ll,i,j,k,iblock)= &
1220 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1229 ! Control printout of the coefficients of virtual-bond-angle potentials
1232 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1237 write (iout,'(//4a)') &
1238 'Type ',onelett(i),onelett(j),onelett(k)
1239 write (iout,'(//a,10x,a)') " l","a[l]"
1240 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1241 write (iout,'(i2,1pe15.5)') &
1242 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1244 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1245 "b",l,"c",l,"d",l,"e",l
1247 write (iout,'(i2,4(1pe15.5))') m,&
1248 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1249 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1253 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1254 "f+",l,"f-",l,"g+",l,"g-",l
1257 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1258 ffthet(n,m,l,i,j,k,iblock),&
1259 ffthet(m,n,l,i,j,k,iblock),&
1260 ggthet(n,m,l,i,j,k,iblock),&
1261 ggthet(m,n,l,i,j,k,iblock)
1272 !C here will be the apropriate recalibrating for D-aminoacid
1273 read (ithep,*,end=121,err=121) nthetyp
1274 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1275 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1276 do i=-nthetyp+1,nthetyp-1
1277 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1278 do j=0,nbend_kcc_Tb(i)
1279 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1283 write (iout,'(a)') &
1284 "Parameters of the valence-only potentials"
1285 do i=-nthetyp+1,nthetyp-1
1286 write (iout,'(2a)') "Type ",toronelet(i)
1287 do k=0,nbend_kcc_Tb(i)
1288 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1294 write (2,*) "Start reading THETA_PDB",ithep_pdb
1296 ! write (2,*) 'i=',i
1297 read (ithep_pdb,*,err=111,end=111) &
1298 a0thet(i),(athet(j,i,1,1),j=1,2),&
1299 (bthet(j,i,1,1),j=1,2)
1300 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1301 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1302 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1303 sigc0(i)=sigc0(i)**2
1306 athet(1,i,1,-1)=athet(1,i,1,1)
1307 athet(2,i,1,-1)=athet(2,i,1,1)
1308 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1309 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1310 athet(1,i,-1,1)=-athet(1,i,1,1)
1311 athet(2,i,-1,1)=-athet(2,i,1,1)
1312 bthet(1,i,-1,1)=bthet(1,i,1,1)
1313 bthet(2,i,-1,1)=bthet(2,i,1,1)
1316 a0thet(i)=a0thet(-i)
1317 athet(1,i,-1,-1)=athet(1,-i,1,1)
1318 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1319 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1320 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1321 athet(1,i,-1,1)=athet(1,-i,1,1)
1322 athet(2,i,-1,1)=-athet(2,-i,1,1)
1323 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1324 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1325 athet(1,i,1,-1)=-athet(1,-i,1,1)
1326 athet(2,i,1,-1)=athet(2,-i,1,1)
1327 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1328 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1329 theta0(i)=theta0(-i)
1333 polthet(j,i)=polthet(j,-i)
1336 gthet(j,i)=gthet(j,-i)
1339 write (2,*) "End reading THETA_PDB"
1343 !--------------- Reading theta parameters for nucleic acid-------
1344 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1345 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1346 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1347 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1348 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1349 nthetyp_nucl+1,nthetyp_nucl+1))
1350 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1351 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1352 nthetyp_nucl+1,nthetyp_nucl+1))
1353 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1354 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1355 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1356 nthetyp_nucl+1,nthetyp_nucl+1))
1357 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1358 nthetyp_nucl+1,nthetyp_nucl+1))
1359 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1360 nthetyp_nucl+1,nthetyp_nucl+1))
1361 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1362 nthetyp_nucl+1,nthetyp_nucl+1))
1363 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1364 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1365 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1366 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1367 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1368 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1370 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1371 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1373 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1375 aa0thet_nucl(:,:,:)=0.0d0
1376 aathet_nucl(:,:,:,:)=0.0d0
1377 bbthet_nucl(:,:,:,:,:)=0.0d0
1378 ccthet_nucl(:,:,:,:,:)=0.0d0
1379 ddthet_nucl(:,:,:,:,:)=0.0d0
1380 eethet_nucl(:,:,:,:,:)=0.0d0
1381 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1382 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1387 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1388 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1389 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1390 read (ithep_nucl,*,end=111,err=111) &
1391 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1392 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1393 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1394 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1395 read (ithep_nucl,*,end=111,err=111) &
1396 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1397 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1398 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1403 !-------------------------------------------
1404 allocate(nlob(ntyp1)) !(ntyp1)
1405 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1406 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1407 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1414 gaussc(:,:,:,:)=0.0D0
1418 ! Read the parameters of the probability distribution/energy expression
1419 ! of the side chains.
1422 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1426 dsc_inv(i)=1.0D0/dsc(i)
1437 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1438 ((blower(k,l,1),l=1,k),k=1,3)
1439 censc(1,1,-i)=censc(1,1,i)
1440 censc(2,1,-i)=censc(2,1,i)
1441 censc(3,1,-i)=-censc(3,1,i)
1443 read (irotam,*,end=112,err=112) bsc(j,i)
1444 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1445 ((blower(k,l,j),l=1,k),k=1,3)
1446 censc(1,j,-i)=censc(1,j,i)
1447 censc(2,j,-i)=censc(2,j,i)
1448 censc(3,j,-i)=-censc(3,j,i)
1449 ! BSC is amplitude of Gaussian
1456 akl=akl+blower(k,m,j)*blower(l,m,j)
1460 if (((k.eq.3).and.(l.ne.3)) &
1461 .or.((l.eq.3).and.(k.ne.3))) then
1462 gaussc(k,l,j,-i)=-akl
1463 gaussc(l,k,j,-i)=-akl
1465 gaussc(k,l,j,-i)=akl
1466 gaussc(l,k,j,-i)=akl
1475 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1478 if (nlobi.gt.0) then
1480 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1481 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1482 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1483 'log h',(bsc(j,i),j=1,nlobi)
1484 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1485 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1487 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1488 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1491 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1492 write (iout,'(a,f10.4,4(16x,f10.4))') &
1493 'Center ',(bsc(j,i),j=1,nlobi)
1494 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1503 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1504 ! added by Urszula Kozlowska 07/11/2007
1506 !el Maximum number of SC local term fitting function coefficiants
1507 !el integer,parameter :: maxsccoef=65
1509 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1512 read (irotam,*,end=112,err=112)
1514 read (irotam,*,end=112,err=112)
1517 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1522 allocate(nterm_scend(2,ntyp))
1523 allocate(arotam_end(0:6,2,ntyp))
1526 read (irotam_end,*) ijunk
1527 !c write (iout,*) "ijunk",ijunk
1531 read (irotam_end,'(a)')
1532 read (irotam_end,*) nterm_scend(j,i)
1533 !c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
1534 do k=0,nterm_scend(j,i)
1535 read (irotam_end,*) ijunk,arotam_end(k,j,i)
1536 !c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
1542 write (iout,'(a)') &
1543 "Parameters of the local potentials of sidechain ends"
1545 write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
1546 restyp(i,1),restyp(i,1)
1547 do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
1548 write (iout,'(i5,2f20.10)') &
1549 j,arotam_end(j,1,i),arotam_end(j,2,i)
1556 !---------reading nucleic acid parameters for rotamers-------------------
1557 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1558 do i=1,ntyp_molec(2)
1559 read (irotam_nucl,*,end=112,err=112)
1561 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1567 write (iout,*) "Base rotamer parameters"
1568 do i=1,ntyp_molec(2)
1569 write (iout,'(a)') restyp(i,2)
1570 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1575 ! Read the parameters of the probability distribution/energy expression
1576 ! of the side chains.
1578 write (2,*) "Start reading ROTAM_PDB"
1580 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1584 dsc_inv(i)=1.0D0/dsc(i)
1595 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1596 ((blower(k,l,1),l=1,k),k=1,3)
1598 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1599 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1600 ((blower(k,l,j),l=1,k),k=1,3)
1607 akl=akl+blower(k,m,j)*blower(l,m,j)
1617 write (2,*) "End reading ROTAM_PDB"
1623 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1624 !C interaction energy of the Gly, Ala, and Pro prototypes.
1626 read (ifourier,*) nloctyp
1627 SPLIT_FOURIERTOR = nloctyp.lt.0
1628 nloctyp = iabs(nloctyp)
1629 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1630 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1631 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1632 !C allocate(ctilde(2,2,nres))
1633 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1634 !C allocate(gtb1(2,nres))
1635 !C allocate(gtb2(2,nres))
1636 !C allocate(cc(2,2,nres))
1637 !C allocate(dd(2,2,nres))
1638 !C allocate(ee(2,2,nres))
1639 !C allocate(gtcc(2,2,nres))
1640 !C allocate(gtdd(2,2,nres))
1641 !C allocate(gtee(2,2,nres))
1644 allocate(itype2loc(-ntyp1:ntyp1))
1645 allocate(iloctyp(-nloctyp:nloctyp))
1646 allocate(bnew1(3,2,-nloctyp:nloctyp))
1647 allocate(bnew2(3,2,-nloctyp:nloctyp))
1648 allocate(ccnew(3,2,-nloctyp:nloctyp))
1649 allocate(ddnew(3,2,-nloctyp:nloctyp))
1650 allocate(e0new(3,-nloctyp:nloctyp))
1651 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1652 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1653 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1654 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1655 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1656 allocate(e0newtor(3,-nloctyp:nloctyp))
1657 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1670 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1671 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1672 itype2loc(ntyp1)=nloctyp
1673 iloctyp(nloctyp)=ntyp1
1675 itype2loc(-i)=-itype2loc(i)
1678 allocate(iloctyp(-nloctyp:nloctyp))
1679 allocate(itype2loc(-ntyp1:ntyp1))
1686 iloctyp(-i)=-iloctyp(i)
1688 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1689 !c write (iout,*) "nloctyp",nloctyp,
1690 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1691 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1692 !c write (iout,*) "nloctyp",nloctyp,
1693 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1696 !c write (iout,*) "NEWCORR",i
1697 read (ifourier,*,end=115,err=115)
1700 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1703 !c write (iout,*) "NEWCORR BNEW1"
1704 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1707 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1710 !c write (iout,*) "NEWCORR BNEW2"
1711 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1713 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1714 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1716 !c write (iout,*) "NEWCORR CCNEW"
1717 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1719 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1720 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1722 !c write (iout,*) "NEWCORR DDNEW"
1723 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1727 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1731 !c write (iout,*) "NEWCORR EENEW1"
1732 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1734 read (ifourier,*,end=115,err=115) e0new(ii,i)
1736 !c write (iout,*) (e0new(ii,i),ii=1,3)
1738 !c write (iout,*) "NEWCORR EENEW"
1741 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1742 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1743 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1744 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1749 bnew1(ii,1,-i)= bnew1(ii,1,i)
1750 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1751 bnew2(ii,1,-i)= bnew2(ii,1,i)
1752 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1755 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1756 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1757 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1758 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1759 ccnew(ii,1,-i)=ccnew(ii,1,i)
1760 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1761 ddnew(ii,1,-i)=ddnew(ii,1,i)
1762 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1764 e0new(1,-i)= e0new(1,i)
1765 e0new(2,-i)=-e0new(2,i)
1766 e0new(3,-i)=-e0new(3,i)
1768 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1769 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1770 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1771 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1775 write (iout,'(a)') "Coefficients of the multibody terms"
1776 do i=-nloctyp+1,nloctyp-1
1777 write (iout,*) "Type: ",onelet(iloctyp(i))
1778 write (iout,*) "Coefficients of the expansion of B1"
1780 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1782 write (iout,*) "Coefficients of the expansion of B2"
1784 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1786 write (iout,*) "Coefficients of the expansion of C"
1787 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1788 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1789 write (iout,*) "Coefficients of the expansion of D"
1790 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1791 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1792 write (iout,*) "Coefficients of the expansion of E"
1793 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1796 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1801 IF (SPLIT_FOURIERTOR) THEN
1803 !c write (iout,*) "NEWCORR TOR",i
1804 read (ifourier,*,end=115,err=115)
1807 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1810 !c write (iout,*) "NEWCORR BNEW1 TOR"
1811 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1814 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1817 !c write (iout,*) "NEWCORR BNEW2 TOR"
1818 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1820 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1821 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1823 !c write (iout,*) "NEWCORR CCNEW TOR"
1824 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1826 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1827 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1829 !c write (iout,*) "NEWCORR DDNEW TOR"
1830 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1834 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1838 !c write (iout,*) "NEWCORR EENEW1 TOR"
1839 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1841 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1843 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1845 !c write (iout,*) "NEWCORR EENEW TOR"
1848 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1849 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1850 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1851 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1856 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1857 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1858 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1859 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1862 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1863 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1864 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1865 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1867 e0newtor(1,-i)= e0newtor(1,i)
1868 e0newtor(2,-i)=-e0newtor(2,i)
1869 e0newtor(3,-i)=-e0newtor(3,i)
1871 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1872 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1873 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1874 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1878 write (iout,'(a)') &
1879 "Single-body coefficients of the torsional potentials"
1880 do i=-nloctyp+1,nloctyp-1
1881 write (iout,*) "Type: ",onelet(iloctyp(i))
1882 write (iout,*) "Coefficients of the expansion of B1tor"
1884 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1885 j,(bnew1tor(k,j,i),k=1,3)
1887 write (iout,*) "Coefficients of the expansion of B2tor"
1889 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1890 j,(bnew2tor(k,j,i),k=1,3)
1892 write (iout,*) "Coefficients of the expansion of Ctor"
1893 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1894 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1895 write (iout,*) "Coefficients of the expansion of Dtor"
1896 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1897 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1898 write (iout,*) "Coefficients of the expansion of Etor"
1899 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1902 write (iout,'(1hE,2i1,2f10.5)') &
1903 j,k,(eenewtor(l,j,k,i),l=1,2)
1909 do i=-nloctyp+1,nloctyp-1
1912 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1917 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1921 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1922 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1923 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1924 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1927 ENDIF !SPLIT_FOURIER_TOR
1929 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1930 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1931 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1932 allocate(b(13,-nloctyp-1:nloctyp+1))
1934 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1936 read (ifourier,*,end=115,err=115)
1937 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1939 write (iout,*) 'Type ',onelet(iloctyp(i))
1940 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1954 !c B1tilde(1,i) = b(3)
1955 !c! B1tilde(2,i) =-b(5)
1956 !c! B1tilde(1,-i) =-b(3)
1957 !c! B1tilde(2,-i) =b(5)
1958 !c! b1tilde(1,i)=0.0d0
1959 !c b1tilde(2,i)=0.0d0
1964 !cc B1tilde(1,i) = b(3,i)
1965 !cc B1tilde(2,i) =-b(5,i)
1966 !c B1tilde(1,-i) =-b(3,i)
1967 !c B1tilde(2,-i) =b(5,i)
1968 !cc b1tilde(1,i)=0.0d0
1969 !cc b1tilde(2,i)=0.0d0
1970 !cc B2(1,i) = b(2,i)
1971 !cc B2(2,i) = b(4,i)
1973 !c B2(2,-i) =-b(4,i)
1977 CCold(1,1,i)= b(7,i)
1978 CCold(2,2,i)=-b(7,i)
1979 CCold(2,1,i)= b(9,i)
1980 CCold(1,2,i)= b(9,i)
1981 CCold(1,1,-i)= b(7,i)
1982 CCold(2,2,-i)=-b(7,i)
1983 CCold(2,1,-i)=-b(9,i)
1984 CCold(1,2,-i)=-b(9,i)
1989 !c Ctilde(1,1,i)= CCold(1,1,i)
1990 !c Ctilde(1,2,i)= CCold(1,2,i)
1991 !c Ctilde(2,1,i)=-CCold(2,1,i)
1992 !c Ctilde(2,2,i)=-CCold(2,2,i)
1997 !c Ctilde(1,1,i)= CCold(1,1,i)
1998 !c Ctilde(1,2,i)= CCold(1,2,i)
1999 !c Ctilde(2,1,i)=-CCold(2,1,i)
2000 !c Ctilde(2,2,i)=-CCold(2,2,i)
2002 !c Ctilde(1,1,i)=0.0d0
2003 !c Ctilde(1,2,i)=0.0d0
2004 !c Ctilde(2,1,i)=0.0d0
2005 !c Ctilde(2,2,i)=0.0d0
2006 DDold(1,1,i)= b(6,i)
2007 DDold(2,2,i)=-b(6,i)
2008 DDold(2,1,i)= b(8,i)
2009 DDold(1,2,i)= b(8,i)
2010 DDold(1,1,-i)= b(6,i)
2011 DDold(2,2,-i)=-b(6,i)
2012 DDold(2,1,-i)=-b(8,i)
2013 DDold(1,2,-i)=-b(8,i)
2018 !c Dtilde(1,1,i)= DD(1,1,i)
2019 !c Dtilde(1,2,i)= DD(1,2,i)
2020 !c Dtilde(2,1,i)=-DD(2,1,i)
2021 !c Dtilde(2,2,i)=-DD(2,2,i)
2023 !c Dtilde(1,1,i)=0.0d0
2024 !c Dtilde(1,2,i)=0.0d0
2025 !c Dtilde(2,1,i)=0.0d0
2026 !c Dtilde(2,2,i)=0.0d0
2027 EEold(1,1,i)= b(10,i)+b(11,i)
2028 EEold(2,2,i)=-b(10,i)+b(11,i)
2029 EEold(2,1,i)= b(12,i)-b(13,i)
2030 EEold(1,2,i)= b(12,i)+b(13,i)
2031 EEold(1,1,-i)= b(10,i)+b(11,i)
2032 EEold(2,2,-i)=-b(10,i)+b(11,i)
2033 EEold(2,1,-i)=-b(12,i)+b(13,i)
2034 EEold(1,2,-i)=-b(12,i)-b(13,i)
2035 write(iout,*) "TU DOCHODZE"
2041 !c ee(2,1,i)=ee(1,2,i)
2046 "Coefficients of the cumulants (independent of valence angles)"
2047 do i=-nloctyp+1,nloctyp-1
2048 write (iout,*) 'Type ',onelet(iloctyp(i))
2050 write(iout,'(2f10.5)') B(3,i),B(5,i)
2052 write(iout,'(2f10.5)') B(2,i),B(4,i)
2055 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
2059 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
2063 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
2072 ! Read torsional parameters in old format
2074 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
2076 read (itorp,*,end=113,err=113) ntortyp,nterm_old
2077 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
2078 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2080 !el from energy module--------
2081 allocate(v1(nterm_old,ntortyp,ntortyp))
2082 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2083 !el---------------------------
2088 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2094 write (iout,'(/a/)') 'Torsional constants:'
2097 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2098 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2104 ! Read torsional parameters
2106 IF (TOR_MODE.eq.0) THEN
2107 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2108 read (itorp,*,end=113,err=113) ntortyp
2109 !el from energy module---------
2110 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2113 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2114 allocate(vlor2(maxlor,ntortyp,ntortyp))
2115 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2116 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2118 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2119 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2120 !el---------------------------
2123 !el---------------------------
2125 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2127 itortyp(i)=-itortyp(-i)
2129 itortyp(ntyp1)=ntortyp
2130 itortyp(-ntyp1)=-ntortyp
2132 write (iout,*) 'ntortyp',ntortyp
2134 do j=-ntortyp+1,ntortyp-1
2135 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2137 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2138 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2141 do k=1,nterm(i,j,iblock)
2142 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2144 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2145 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2146 v0ij=v0ij+si*v1(k,i,j,iblock)
2148 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2149 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2150 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2152 do k=1,nlor(i,j,iblock)
2153 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2154 vlor2(k,i,j),vlor3(k,i,j)
2155 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2158 v0(-i,-j,iblock)=v0ij
2164 write (iout,'(/a/)') 'Torsional constants:'
2166 do i=-ntortyp,ntortyp
2167 do j=-ntortyp,ntortyp
2168 write (iout,*) 'ityp',i,' jtyp',j
2169 write (iout,*) 'Fourier constants'
2170 do k=1,nterm(i,j,iblock)
2171 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2174 write (iout,*) 'Lorenz constants'
2175 do k=1,nlor(i,j,iblock)
2176 write (iout,'(3(1pe15.5))') &
2177 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2183 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2185 ! 6/23/01 Read parameters for double torsionals
2187 !el from energy module------------
2188 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2189 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2190 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2191 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2192 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2193 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2194 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2195 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2196 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2197 !---------------------------------
2201 do j=-ntortyp+1,ntortyp-1
2202 do k=-ntortyp+1,ntortyp-1
2203 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2204 ! write (iout,*) "OK onelett",
2207 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2208 .or. t3.ne.toronelet(k)) then
2209 write (iout,*) "Error in double torsional parameter file",&
2212 call MPI_Finalize(Ierror)
2214 stop "Error in double torsional parameter file"
2216 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2217 ntermd_2(i,j,k,iblock)
2218 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2219 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2220 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2221 ntermd_1(i,j,k,iblock))
2222 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2223 ntermd_1(i,j,k,iblock))
2224 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2225 ntermd_1(i,j,k,iblock))
2226 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2227 ntermd_1(i,j,k,iblock))
2228 ! Martix of D parameters for one dimesional foureir series
2229 do l=1,ntermd_1(i,j,k,iblock)
2230 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2231 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2232 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2233 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2234 ! write(iout,*) "whcodze" ,
2235 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2237 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2238 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2239 v2s(m,l,i,j,k,iblock),&
2240 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2241 ! Martix of D parameters for two dimesional fourier series
2242 do l=1,ntermd_2(i,j,k,iblock)
2244 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2245 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2246 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2247 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2256 write (iout,*) 'Constants for double torsionals'
2259 do j=-ntortyp+1,ntortyp-1
2260 do k=-ntortyp+1,ntortyp-1
2261 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2262 ' nsingle',ntermd_1(i,j,k,iblock),&
2263 ' ndouble',ntermd_2(i,j,k,iblock)
2265 write (iout,*) 'Single angles:'
2266 do l=1,ntermd_1(i,j,k,iblock)
2267 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2268 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2269 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2270 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2273 write (iout,*) 'Pairs of angles:'
2274 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2275 do l=1,ntermd_2(i,j,k,iblock)
2276 write (iout,'(i5,20f10.5)') &
2277 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2280 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2281 do l=1,ntermd_2(i,j,k,iblock)
2282 write (iout,'(i5,20f10.5)') &
2283 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2284 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2294 itype2loc(i)=itortyp(i)
2298 ELSE IF (TOR_MODE.eq.1) THEN
2300 !C read valence-torsional parameters
2301 read (itorp,*,end=121,err=121) ntortyp
2303 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2305 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2306 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2309 itype2loc(i)=itortyp(i)
2313 itortyp(i)=-itortyp(-i)
2315 do i=-ntortyp+1,ntortyp-1
2316 do j=-ntortyp+1,ntortyp-1
2317 !C first we read the cos and sin gamma parameters
2318 read (itorp,'(13x,a)',end=121,err=121) string
2319 write (iout,*) i,j,string
2320 read (itorp,*,end=121,err=121) &
2321 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2322 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2323 do k=1,nterm_kcc(j,i)
2324 do l=1,nterm_kcc_Tb(j,i)
2325 do ll=1,nterm_kcc_Tb(j,i)
2326 read (itorp,*,end=121,err=121) ii,jj,kk, &
2327 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2335 !c AL 4/8/16: Calculate coefficients from one-body parameters
2337 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2338 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2339 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2340 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2341 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2344 print *,i,itortyp(i)
2345 itortyp(i)=itype2loc(i)
2348 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2350 do i=-ntortyp+1,ntortyp-1
2351 do j=-ntortyp+1,ntortyp-1
2354 do k=1,nterm_kcc_Tb(j,i)
2355 do l=1,nterm_kcc_Tb(j,i)
2356 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2357 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2358 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2359 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2362 do k=1,nterm_kcc_Tb(j,i)
2363 do l=1,nterm_kcc_Tb(j,i)
2365 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2366 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2367 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2368 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2370 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2371 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2372 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2373 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2377 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2381 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2386 if (tor_mode.gt.0 .and. lprint) then
2387 !c Print valence-torsional parameters
2388 write (iout,'(a)') &
2389 "Parameters of the valence-torsional potentials"
2390 do i=-ntortyp+1,ntortyp-1
2391 do j=-ntortyp+1,ntortyp-1
2392 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2393 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2394 do k=1,nterm_kcc(j,i)
2395 do l=1,nterm_kcc_Tb(j,i)
2396 do ll=1,nterm_kcc_Tb(j,i)
2397 write (iout,'(3i5,2f15.4)')&
2398 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2406 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2407 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2408 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2409 !el from energy module---------
2410 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2411 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2413 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2414 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2415 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2416 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2418 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2419 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2420 !el---------------------------
2423 !el--------------------
2424 read (itorp_nucl,*,end=113,err=113) &
2425 (itortyp_nucl(i),i=1,ntyp_molec(2))
2426 ! print *,itortyp_nucl(:)
2427 !c write (iout,*) 'ntortyp',ntortyp
2430 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2431 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2434 do k=1,nterm_nucl(i,j)
2435 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2436 v0ij=v0ij+si*v1_nucl(k,i,j)
2439 do k=1,nlor_nucl(i,j)
2440 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2441 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2442 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2448 ! Read of Side-chain backbone correlation parameters
2449 ! Modified 11 May 2012 by Adasko
2452 read (isccor,*,end=119,err=119) nsccortyp
2454 !el from module energy-------------
2455 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2456 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2457 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2458 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2459 !-----------------------------------
2461 !el from module energy-------------
2462 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2464 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2466 isccortyp(i)=-isccortyp(-i)
2468 iscprol=isccortyp(20)
2469 ! write (iout,*) 'ntortyp',ntortyp
2471 !c maxinter is maximum interaction sites
2472 !el from module energy---------
2473 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2474 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2475 -nsccortyp:nsccortyp))
2476 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2477 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2478 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2479 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2480 !-----------------------------------
2484 read (isccor,*,end=119,err=119) &
2485 nterm_sccor(i,j),nlor_sccor(i,j)
2491 nterm_sccor(-i,j)=nterm_sccor(i,j)
2492 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2493 nterm_sccor(i,-j)=nterm_sccor(i,j)
2494 do k=1,nterm_sccor(i,j)
2495 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2497 if (j.eq.iscprol) then
2498 if (i.eq.isccortyp(10)) then
2499 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2500 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2502 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2503 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2504 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2505 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2506 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2507 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2508 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2509 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2512 if (i.eq.isccortyp(10)) then
2513 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2514 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2516 if (j.eq.isccortyp(10)) then
2517 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2518 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2520 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2521 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2522 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2523 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2524 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2525 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2529 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2530 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2531 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2532 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2535 do k=1,nlor_sccor(i,j)
2536 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2537 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2538 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2539 (1+vlor3sccor(k,i,j)**2)
2541 v0sccor(l,i,j)=v0ijsccor
2542 v0sccor(l,-i,j)=v0ijsccor1
2543 v0sccor(l,i,-j)=v0ijsccor2
2544 v0sccor(l,-i,-j)=v0ijsccor3
2550 !el from module energy-------------
2551 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2553 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2554 ! write (iout,*) 'ntortyp',ntortyp
2556 !c maxinter is maximum interaction sites
2557 !el from module energy---------
2558 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2559 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2560 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2561 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2562 !-----------------------------------
2566 read (isccor,*,end=119,err=119) &
2567 nterm_sccor(i,j),nlor_sccor(i,j)
2571 do k=1,nterm_sccor(i,j)
2572 print *,"test",i,j,k,l
2573 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2575 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2578 do k=1,nlor_sccor(i,j)
2579 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2580 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2581 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2582 (1+vlor3sccor(k,i,j)**2)
2584 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2593 write (iout,'(/a/)') 'Torsional constants:'
2596 write (iout,*) 'ityp',i,' jtyp',j
2597 write (iout,*) 'Fourier constants'
2598 do k=1,nterm_sccor(i,j)
2599 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2601 write (iout,*) 'Lorenz constants'
2602 do k=1,nlor_sccor(i,j)
2603 write (iout,'(3(1pe15.5))') &
2604 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2611 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2612 ! interaction energy of the Gly, Ala, and Pro prototypes.
2615 ! Read electrostatic-interaction parameters
2620 write (iout,'(/a)') 'Electrostatic interaction constants:'
2621 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2622 'IT','JT','APP','BPP','AEL6','AEL3'
2624 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2625 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2626 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2627 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2632 app (i,j)=epp(i,j)*rri*rri
2633 bpp (i,j)=-2.0D0*epp(i,j)*rri
2634 ael6(i,j)=elpp6(i,j)*4.2D0**6
2635 ael3(i,j)=elpp3(i,j)*4.2D0**3
2637 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2643 ! Read side-chain interaction parameters.
2645 !el from module energy - COMMON.INTERACT-------
2646 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2647 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2648 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2650 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2651 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2652 allocate(epslip(ntyp,ntyp))
2660 !--------------------------------
2662 read (isidep,*,end=117,err=117) ipot,expon
2663 if (ipot.lt.1 .or. ipot.gt.5) then
2664 write (iout,'(2a)') 'Error while reading SC interaction',&
2665 'potential file - unknown potential type.'
2667 call MPI_Finalize(Ierror)
2673 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2674 ', exponents are ',expon,2*expon
2675 ! goto (10,20,30,30,40) ipot
2677 !----------------------- LJ potential ---------------------------------
2679 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2680 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2681 (sigma0(i),i=1,ntyp)
2683 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2684 write (iout,'(a/)') 'The epsilon array:'
2685 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2686 write (iout,'(/a)') 'One-body parameters:'
2687 write (iout,'(a,4x,a)') 'residue','sigma'
2688 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2691 !----------------------- LJK potential --------------------------------
2693 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2694 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2695 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2697 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2698 write (iout,'(a/)') 'The epsilon array:'
2699 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2700 write (iout,'(/a)') 'One-body parameters:'
2701 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2702 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2706 !---------------------- GB or BP potential -----------------------------
2709 ! print *,"I AM in SCELE",scelemode
2710 if (scelemode.eq.0) then
2712 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2714 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2715 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2716 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2717 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2719 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2722 ! For the GB potential convert sigma'**2 into chi'
2725 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2729 write (iout,'(/a/)') 'Parameters of the BP potential:'
2730 write (iout,'(a/)') 'The epsilon array:'
2731 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2732 write (iout,'(/a)') 'One-body parameters:'
2733 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2735 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2736 chip(i),alp(i),i=1,ntyp)
2739 ! print *,ntyp,"NTYP"
2740 allocate(icharge(ntyp1))
2741 ! print *,ntyp,icharge(i)
2743 read (isidep,*) (icharge(i),i=1,ntyp)
2744 print *,ntyp,icharge(i)
2745 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2746 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2747 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2748 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2749 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2750 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2751 allocate(epsintab(ntyp,ntyp))
2752 allocate(dtail(2,ntyp,ntyp))
2753 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2754 allocate(wqdip(2,ntyp,ntyp))
2755 allocate(wstate(4,ntyp,ntyp))
2756 allocate(dhead(2,2,ntyp,ntyp))
2757 allocate(nstate(ntyp,ntyp))
2758 allocate(debaykap(ntyp,ntyp))
2760 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2761 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2765 ! write (*,*) "Im in ALAB", i, " ", j
2767 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2768 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2769 chis(i,j),chis(j,i), & !2 w tej linii
2770 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2771 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2772 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2773 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2774 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2775 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2776 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2777 IF ((LaTeX).and.(i.gt.24)) then
2778 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2779 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2780 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2781 chis(i,j),chis(j,i) !2 w tej linii
2783 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2787 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2791 IF ((LaTeX).and.(i.gt.24)) then
2792 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2793 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2794 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2795 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2796 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2797 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2798 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2805 sigma(i,j) = sigma(j,i)
2806 sigmap1(i,j)=sigmap1(j,i)
2807 sigmap2(i,j)=sigmap2(j,i)
2808 sigiso1(i,j)=sigiso1(j,i)
2809 sigiso2(i,j)=sigiso2(j,i)
2810 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2811 nstate(i,j) = nstate(j,i)
2812 dtail(1,i,j) = dtail(2,j,i)
2813 dtail(2,i,j) = dtail(1,j,i)
2815 alphasur(k,i,j) = alphasur(k,j,i)
2816 wstate(k,i,j) = wstate(k,j,i)
2817 alphiso(k,i,j) = alphiso(k,j,i)
2820 dhead(2,1,i,j) = dhead(1,1,j,i)
2821 dhead(2,2,i,j) = dhead(1,2,j,i)
2822 dhead(1,1,i,j) = dhead(2,1,j,i)
2823 dhead(1,2,i,j) = dhead(2,2,j,i)
2825 epshead(i,j) = epshead(j,i)
2826 sig0head(i,j) = sig0head(j,i)
2829 wqdip(k,i,j) = wqdip(k,j,i)
2832 wquad(i,j) = wquad(j,i)
2833 epsintab(i,j) = epsintab(j,i)
2834 debaykap(i,j)=debaykap(j,i)
2835 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2842 !--------------------- GBV potential -----------------------------------
2844 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2845 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2846 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2847 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2849 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2850 write (iout,'(a/)') 'The epsilon array:'
2851 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2852 write (iout,'(/a)') 'One-body parameters:'
2853 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2854 's||/s_|_^2',' chip ',' alph '
2855 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2856 sigii(i),chip(i),alp(i),i=1,ntyp)
2859 write(iout,*)"Wrong ipot"
2865 !-----------------------------------------------------------------------
2866 ! Calculate the "working" parameters of SC interactions.
2868 !el from module energy - COMMON.INTERACT-------
2869 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2870 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2871 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2872 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2873 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2874 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2876 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2882 if (scelemode.eq.0) then
2891 sc_aa_tube_par(:)=0.0d0
2892 sc_bb_tube_par(:)=0.0d0
2894 !--------------------------------
2899 epslip(i,j)=epslip(j,i)
2902 if (scelemode.eq.0) then
2905 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2906 sigma(j,i)=sigma(i,j)
2907 rs0(i,j)=dwa16*sigma(i,j)
2912 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2913 'Working parameters of the SC interactions:',&
2914 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2919 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2921 ! print *,"SIGMA ZLA?",sigma(i,j)
2929 sigeps=dsign(1.0D0,epsij)
2931 aa_aq(i,j)=epsij*rrij*rrij
2932 ! print *,"ADASKO",epsij,rrij,expon
2933 bb_aq(i,j)=-sigeps*epsij*rrij
2934 aa_aq(j,i)=aa_aq(i,j)
2935 bb_aq(j,i)=bb_aq(i,j)
2936 epsijlip=epslip(i,j)
2937 sigeps=dsign(1.0D0,epsijlip)
2938 epsijlip=dabs(epsijlip)
2939 aa_lip(i,j)=epsijlip*rrij*rrij
2940 bb_lip(i,j)=-sigeps*epsijlip*rrij
2941 aa_lip(j,i)=aa_lip(i,j)
2942 bb_lip(j,i)=bb_lip(i,j)
2943 !C write(iout,*) aa_lip
2944 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2945 sigt1sq=sigma0(i)**2
2946 sigt2sq=sigma0(j)**2
2949 ratsig1=sigt2sq/sigt1sq
2950 ratsig2=1.0D0/ratsig1
2951 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2952 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2953 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2957 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2958 sigmaii(i,j)=rsum_max
2959 sigmaii(j,i)=rsum_max
2961 ! sigmaii(i,j)=r0(i,j)
2962 ! sigmaii(j,i)=r0(i,j)
2964 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2965 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2966 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2967 augm(i,j)=epsij*r_augm**(2*expon)
2968 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2975 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2976 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2977 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2982 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2983 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2984 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2985 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2986 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2987 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2988 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2989 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2990 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2991 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2992 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2993 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2994 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2995 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2996 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2997 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
3006 read (isidep_nucl,*) ipot_nucl
3007 ! print *,"TU?!",ipot_nucl
3008 if (ipot_nucl.eq.1) then
3009 do i=1,ntyp_molec(2)
3010 do j=i,ntyp_molec(2)
3011 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
3012 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
3016 do i=1,ntyp_molec(2)
3017 do j=i,ntyp_molec(2)
3018 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
3019 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
3020 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
3024 ! rpp(1,1)=2**(1.0/6.0)*5.16158
3025 do i=1,ntyp_molec(2)
3026 do j=i,ntyp_molec(2)
3027 rrij=sigma_nucl(i,j)
3031 epsij=4*eps_nucl(i,j)
3032 sigeps=dsign(1.0D0,epsij)
3034 aa_nucl(i,j)=epsij*rrij*rrij
3035 bb_nucl(i,j)=-sigeps*epsij*rrij
3036 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
3037 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
3038 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
3039 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
3040 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
3041 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
3044 aa_nucl(i,j)=aa_nucl(j,i)
3045 bb_nucl(i,j)=bb_nucl(j,i)
3046 ael3_nucl(i,j)=ael3_nucl(j,i)
3047 ael6_nucl(i,j)=ael6_nucl(j,i)
3048 ael63_nucl(i,j)=ael63_nucl(j,i)
3049 ael32_nucl(i,j)=ael32_nucl(j,i)
3050 elpp3_nucl(i,j)=elpp3_nucl(j,i)
3051 elpp6_nucl(i,j)=elpp6_nucl(j,i)
3052 elpp63_nucl(i,j)=elpp63_nucl(j,i)
3053 elpp32_nucl(i,j)=elpp32_nucl(j,i)
3054 eps_nucl(i,j)=eps_nucl(j,i)
3055 sigma_nucl(i,j)=sigma_nucl(j,i)
3056 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
3060 write(iout,*) "tube param"
3061 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
3062 ccavtubpep,dcavtubpep,tubetranenepep
3063 sigmapeptube=sigmapeptube**6
3064 sigeps=dsign(1.0D0,epspeptube)
3065 epspeptube=dabs(epspeptube)
3066 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
3067 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
3068 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
3070 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
3071 ccavtub(i),dcavtub(i),tubetranene(i)
3072 sigmasctube=sigmasctube**6
3073 sigeps=dsign(1.0D0,epssctube)
3074 epssctube=dabs(epssctube)
3075 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
3076 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
3077 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
3079 !-----------------READING SC BASE POTENTIALS-----------------------------
3080 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3081 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3082 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3083 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3084 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3085 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3086 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3087 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3088 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3089 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3090 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3091 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3092 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3093 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3094 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3095 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3097 write (iout,*) "ESCBASEPARM"
3098 do i=1,ntyp_molec(1)
3099 do j=1,ntyp_molec(2) ! without U then we will take T for U
3100 ! write (*,*) "Im in ", i, " ", j
3101 read(isidep_scbase,*) &
3102 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3103 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3104 ! write(*,*) "eps",eps_scbase(i,j)
3105 read(isidep_scbase,*) &
3106 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3107 chis_scbase(i,j,1),chis_scbase(i,j,2)
3108 read(isidep_scbase,*) &
3109 dhead_scbasei(i,j), &
3110 dhead_scbasej(i,j), &
3111 rborn_scbasei(i,j),rborn_scbasej(i,j)
3112 read(isidep_scbase,*) &
3113 (wdipdip_scbase(k,i,j),k=1,3), &
3114 (wqdip_scbase(k,i,j),k=1,2)
3115 read(isidep_scbase,*) &
3116 alphapol_scbase(i,j), &
3117 epsintab_scbase(i,j)
3118 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3119 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3120 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3121 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3122 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3123 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3124 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3125 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3127 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3128 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3129 write(*,*) "eps",eps_scbase(i,j)
3131 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3132 chis_scbase(i,j,1),chis_scbase(i,j,2)
3134 dhead_scbasei(i,j), &
3135 dhead_scbasej(i,j), &
3136 rborn_scbasei(i,j),rborn_scbasej(i,j)
3138 (wdipdip_scbase(k,i,j),k=1,3), &
3139 (wqdip_scbase(k,i,j),k=1,2)
3141 alphapol_scbase(i,j), &
3142 epsintab_scbase(i,j)
3147 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3148 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3149 write(*,*) "eps",eps_scbase(i,j)
3151 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3152 chis_scbase(i,j,1),chis_scbase(i,j,2)
3154 dhead_scbasei(i,j), &
3155 dhead_scbasej(i,j), &
3156 rborn_scbasei(i,j),rborn_scbasej(i,j)
3158 (wdipdip_scbase(k,i,j),k=1,3), &
3159 (wqdip_scbase(k,i,j),k=1,2)
3161 alphapol_scbase(i,j), &
3162 epsintab_scbase(i,j)
3165 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3166 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3168 do i=1,ntyp_molec(1)
3169 do j=1,ntyp_molec(2)-1
3170 epsij=eps_scbase(i,j)
3171 rrij=sigma_scbase(i,j)
3176 sigeps=dsign(1.0D0,epsij)
3178 aa_scbase(i,j)=epsij*rrij*rrij
3179 bb_scbase(i,j)=-sigeps*epsij*rrij
3182 !-----------------READING PEP BASE POTENTIALS-------------------
3183 allocate(eps_pepbase(ntyp_molec(2)))
3184 allocate(sigma_pepbase(ntyp_molec(2)))
3185 allocate(chi_pepbase(ntyp_molec(2),2))
3186 allocate(chipp_pepbase(ntyp_molec(2),2))
3187 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3188 allocate(sigmap1_pepbase(ntyp_molec(2)))
3189 allocate(sigmap2_pepbase(ntyp_molec(2)))
3190 allocate(chis_pepbase(ntyp_molec(2),2))
3191 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3193 write (iout,*) "EPEPBASEPARM"
3195 do j=1,ntyp_molec(2) ! without U then we will take T for U
3196 write (*,*) "Im in ", i, " ", j
3197 read(isidep_pepbase,*) &
3198 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3199 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3200 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3201 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3202 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3203 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3204 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3205 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3206 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3207 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3209 write(*,*) "eps",eps_pepbase(j)
3210 read(isidep_pepbase,*) &
3211 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3212 chis_pepbase(j,1),chis_pepbase(j,2)
3213 read(isidep_pepbase,*) &
3214 (wdipdip_pepbase(k,j),k=1,3)
3215 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3216 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3218 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3219 chis_pepbase(j,1),chis_pepbase(j,2)
3221 (wdipdip_pepbase(k,j),k=1,3)
3225 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3226 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3228 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3229 chis_pepbase(j,1),chis_pepbase(j,2)
3231 (wdipdip_pepbase(k,j),k=1,3)
3233 allocate(aa_pepbase(ntyp_molec(2)))
3234 allocate(bb_pepbase(ntyp_molec(2)))
3236 do j=1,ntyp_molec(2)-1
3237 epsij=eps_pepbase(j)
3238 rrij=sigma_pepbase(j)
3243 sigeps=dsign(1.0D0,epsij)
3245 aa_pepbase(j)=epsij*rrij*rrij
3246 bb_pepbase(j)=-sigeps*epsij*rrij
3248 !--------------READING SC PHOSPHATE-------------------------------------
3249 allocate(eps_scpho(ntyp_molec(1)))
3250 allocate(sigma_scpho(ntyp_molec(1)))
3251 allocate(chi_scpho(ntyp_molec(1),2))
3252 allocate(chipp_scpho(ntyp_molec(1),2))
3253 allocate(alphasur_scpho(4,ntyp_molec(1)))
3254 allocate(sigmap1_scpho(ntyp_molec(1)))
3255 allocate(sigmap2_scpho(ntyp_molec(1)))
3256 allocate(chis_scpho(ntyp_molec(1),2))
3257 allocate(wqq_scpho(ntyp_molec(1)))
3258 allocate(wqdip_scpho(2,ntyp_molec(1)))
3259 allocate(alphapol_scpho(ntyp_molec(1)))
3260 allocate(epsintab_scpho(ntyp_molec(1)))
3261 allocate(dhead_scphoi(ntyp_molec(1)))
3262 allocate(rborn_scphoi(ntyp_molec(1)))
3263 allocate(rborn_scphoj(ntyp_molec(1)))
3264 allocate(alphi_scpho(ntyp_molec(1)))
3268 do j=1,ntyp_molec(1) ! without U then we will take T for U
3269 write (*,*) "Im in scpho ", i, " ", j
3270 read(isidep_scpho,*) &
3271 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3272 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3273 write(*,*) "eps",eps_scpho(j)
3274 read(isidep_scpho,*) &
3275 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3276 chis_scpho(j,1),chis_scpho(j,2)
3277 read(isidep_scpho,*) &
3278 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3279 read(isidep_scpho,*) &
3280 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3282 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3283 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3284 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3285 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3286 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3287 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3288 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3289 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3293 allocate(aa_scpho(ntyp_molec(1)))
3294 allocate(bb_scpho(ntyp_molec(1)))
3296 do j=1,ntyp_molec(1)
3303 sigeps=dsign(1.0D0,epsij)
3305 aa_scpho(j)=epsij*rrij*rrij
3306 bb_scpho(j)=-sigeps*epsij*rrij
3310 read(isidep_peppho,*) &
3311 eps_peppho,sigma_peppho
3312 read(isidep_peppho,*) &
3313 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3314 read(isidep_peppho,*) &
3315 (wqdip_peppho(k),k=1,2)
3323 sigeps=dsign(1.0D0,epsij)
3325 aa_peppho=epsij*rrij*rrij
3326 bb_peppho=-sigeps*epsij*rrij
3329 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3334 ! Define the SC-p interaction constants (hard-coded; old style)
3337 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3339 ! aad(i,1)=0.3D0*4.0D0**12
3340 ! Following line for constants currently implemented
3341 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3342 aad(i,1)=1.5D0*4.0D0**12
3343 ! aad(i,1)=0.17D0*5.6D0**12
3345 ! "Soft" SC-p repulsion
3347 ! Following line for constants currently implemented
3348 ! aad(i,1)=0.3D0*4.0D0**6
3349 ! "Hard" SC-p repulsion
3350 bad(i,1)=3.0D0*4.0D0**6
3351 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3360 ! 8/9/01 Read the SC-p interaction constants from file
3363 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3366 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3367 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3368 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3369 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3373 write (iout,*) "Parameters of SC-p interactions:"
3375 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3376 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3381 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3383 do i=1,ntyp_molec(2)
3384 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3386 do i=1,ntyp_molec(2)
3387 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3388 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3390 r0pp=1.12246204830937298142*5.16158
3396 ! Define the constants of the disulfide bridge
3400 ! Old arbitrary potential - commented out.
3405 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3406 ! energy surface of diethyl disulfide.
3407 ! A. Liwo and U. Kozlowska, 11/24/03
3423 allocate(ichargelipid(ntyp_molec(4)))
3424 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3425 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3426 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3427 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3428 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3429 kjtokcal=0.2390057361
3431 !HERE THE MASS of MARTINI
3432 write(*,*) "before MARTINI PARAM"
3433 do i=1,ntyp_molec(4)
3439 msc(ntyp_molec(4)+1,4)=0.1d0
3440 !relative dielectric constant = 15 for implicit screening
3441 k_coulomb_lip=332.0d0/15.0d0
3442 !kbond = 1250 kJ/(mol*nm*2)
3443 kbondlip=1250.0d0*kjtokcal/100.0d0
3445 lip_angle_force=0.0d0
3446 if (DRY_MARTINI.gt.0) then
3447 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3448 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3449 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3450 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3451 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3452 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3453 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3454 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3456 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3457 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3458 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3459 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3460 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3461 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3462 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3463 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3465 lip_angle_angle=0.0d0
3466 lip_angle_angle(3,12,12)=120.0/krad
3467 lip_angle_angle(3,12,18)=180.0/krad
3468 lip_angle_angle(3,18,16)=180.0/krad
3469 lip_angle_angle(12,18,16)=180.0/krad
3470 lip_angle_angle(18,16,18)=120.0/krad
3471 lip_angle_angle(16,18,18)=180.0/krad
3472 lip_angle_angle(12,18,18)=180.0/krad
3473 lip_angle_angle(18,18,18)=180.0/krad
3474 read(ilipbond,*) temp1
3476 read(ilipbond,*) temp1, lip_bond(i,1), &
3477 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3478 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3479 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3480 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3483 lip_bond(i,j)=lip_bond(i,j)*10
3487 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3488 read(ilipnonbond,*) temp1
3490 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3491 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3492 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3493 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3494 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3496 ! write(*,*) i, lip_eps(i,18)
3498 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3501 read(ilipnonbond,*) temp1
3503 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3504 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3505 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3506 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3507 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3510 lip_sig(i,j)=lip_sig(i,j)*10.0
3513 write(*,*) "after MARTINI PARAM"
3517 allocate(alphapolcat(ntyp,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5)))
3518 allocate(alphapolcat2(ntyp,-1:ntyp_molec(5)))
3519 allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5)))
3520 allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5)))
3521 allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5)))
3522 allocate(epsintabcat(ntyp,-1:ntyp_molec(5)))
3523 allocate(dtailcat(2,ntyp,-1:ntyp_molec(5)))
3524 allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5)))
3525 allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5)))
3526 allocate(wstatecat(4,ntyp,-1:ntyp_molec(5)))
3527 allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5)))
3528 allocate(nstatecat(ntyp,-1:ntyp_molec(5)))
3529 allocate(debaykapcat(ntyp,-1:ntyp_molec(5)))
3531 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1))
3532 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1))
3533 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3534 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
3535 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
3538 if (.not.allocated(ichargecat))&
3539 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
3540 write(*,*) "before ions",oldion
3543 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3544 if (oldion.eq.0) then
3545 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3546 allocate(icharge(1:ntyp1))
3547 read(iion,*) (icharge(i),i=1,ntyp)
3551 print *,ntyp_molec(5)
3552 do i=-ntyp_molec(5),ntyp_molec(5)
3553 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3554 print *,msc(i,5),restok(i,5)
3560 do j=-1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3563 ! do j=1,ntyp_molec(5)
3564 ! write (*,*) "Im in ALAB", i, " ", j
3566 epscat(i,j),sigmacat(i,j), &
3567 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3568 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), & !6
3570 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12
3571 ! chiscat(i,j),chiscat(j,i), &
3572 chis1cat(i,j),chis2cat(i,j), &
3574 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19 !5 w tej lini - 1 integer pierwszy
3575 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&!23
3576 dtailcat(1,i,j),dtailcat(2,i,j), &
3577 epsheadcat(i,j),sig0headcat(i,j), &!27
3579 ! rborncat(i,j),rborncat(j,i),&
3580 rborn1cat(i,j),rborn2cat(i,j),&
3581 (wqdipcat(k,i,j),k=1,2), &!31
3582 alphapolcat(i,j),alphapolcat2(j,i), &!33
3583 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3585 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3586 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3588 ! write (iout,*) 'i= ', i, ' j= ', j
3589 ! write (iout,*) 'epsi0= ', epscat(1,j)
3590 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3591 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3592 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3593 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3594 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3595 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3596 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3597 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3598 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3599 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3600 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3601 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3602 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3603 ! write (iout,*) 'a1= ', rborncat(j,1)
3604 ! write (iout,*) 'a2= ', rborncat(1,j)
3605 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3606 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3607 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3608 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3612 ! If ((i.eq.1).and.(j.eq.27)) then
3613 ! write (iout,*) 'SEP'
3614 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3615 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3620 allocate(aa_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)),&
3621 bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)))
3623 do j=-1,ntyp_molec(5)
3628 sigeps=dsign(1.0D0,epsij)
3630 aa_aq_cat(i,j)=epsij*rrij*rrij
3631 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3636 do j=1,ntyp_molec(5)-1
3638 write (iout,*) 'i= ', i, ' j= ', j
3639 write (iout,*) 'epsi0= ', epscat(i,j)
3640 write (iout,*) 'sigma0= ', sigmacat(i,j)
3641 write (iout,*) 'chi1= ', chi1cat(i,j)
3642 write (iout,*) 'chi1= ', chi2cat(i,j)
3643 write (iout,*) 'chip1= ', chipp1cat(i,j)
3644 write (iout,*) 'chip2= ', chipp2cat(i,j)
3645 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3646 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3647 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3648 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3649 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3650 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3651 write (iout,*) 'chis1= ', chis1cat(i,j)
3652 write (iout,*) 'chis1= ', chis2cat(i,j)
3653 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3654 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3655 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3656 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3657 write (iout,*) 'a1= ', rborn1cat(i,j)
3658 write (iout,*) 'a2= ', rborn2cat(i,j)
3659 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3660 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3661 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3662 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3663 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3664 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3667 If ((i.eq.1).and.(j.eq.27)) then
3668 write (iout,*) 'SEP'
3669 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3670 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3677 ! read number of Zn2+
3678 ! here two denotes the Zn2+ and Cu2+
3679 write(iout,*) "before TRANPARM"
3680 allocate(aomicattr(0:3,2))
3681 allocate(athetacattran(0:6,5,2))
3682 allocate(agamacattran(3,5,2))
3683 allocate(acatshiftdsc(5,2))
3684 allocate(bcatshiftdsc(5,2))
3685 allocate(demorsecat(5,2))
3686 allocate(alphamorsecat(5,2))
3687 allocate(x0catleft(5,2))
3688 allocate(x0catright(5,2))
3689 allocate(x0cattrans(5,2))
3690 allocate(ntrantyp(2))
3691 do i=1,1 ! currently only Zn2+
3693 read(iiontran,*) ntrantyp(i)
3696 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3697 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3698 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3699 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3700 read(iiontran,*) (aomicattr(j,i),j=0,3)
3702 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3703 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3704 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3709 write (iout,'(/a)') "Disulfide bridge parameters:"
3710 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3711 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3712 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3713 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3716 if (shield_mode.gt.0) then
3717 pi=4.0D0*datan(1.0D0)
3718 !C VSolvSphere the volume of solving sphere
3720 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3721 !C there will be no distinction between proline peptide group and normal peptide
3722 !C group in case of shielding parameters
3723 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3724 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3725 write (iout,*) VSolvSphere,VSolvSphere_div
3726 !C long axis of side chain
3728 long_r_sidechain(i)=vbldsc0(1,i)
3729 ! if (scelemode.eq.0) then
3730 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3731 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3733 ! short_r_sidechain(i)=sigma(i,i)
3735 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3742 111 write (iout,*) "Error reading bending energy parameters."
3744 112 write (iout,*) "Error reading rotamer energy parameters."
3746 113 write (iout,*) "Error reading torsional energy parameters."
3748 114 write (iout,*) "Error reading double torsional energy parameters."
3750 115 write (iout,*) &
3751 "Error reading cumulant (multibody energy) parameters."
3753 116 write (iout,*) "Error reading electrostatic energy parameters."
3755 117 write (iout,*) "Error reading side chain interaction parameters."
3757 118 write (iout,*) "Error reading SCp interaction parameters."
3759 119 write (iout,*) "Error reading SCCOR parameters"
3761 121 write (iout,*) "Error in Czybyshev parameters"
3763 123 write(iout,*) "Error in transition metal parameters"
3766 call MPI_Finalize(Ierror)
3770 end subroutine parmread
3772 !-----------------------------------------------------------------------------
3774 !-----------------------------------------------------------------------------
3775 subroutine printmat(ldim,m,n,iout,key,a)
3778 character(len=3),dimension(n) :: key
3779 real(kind=8),dimension(ldim,n) :: a
3781 integer :: i,j,k,m,iout,nlim
3785 write (iout,1000) (key(k),k=i,nlim)
3787 1000 format (/5x,8(6x,a3))
3788 1020 format (/80(1h-)/)
3790 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3793 1010 format (a3,2x,8(f9.4))
3795 end subroutine printmat
3796 !-----------------------------------------------------------------------------
3798 !-----------------------------------------------------------------------------
3800 ! Read the PDB file and convert the peptide geometry into virtual-chain
3803 use energy_data, only: itype,istype
3807 ! use control, only: rescode,sugarcode
3808 ! implicit real*8 (a-h,o-z)
3809 ! include 'DIMENSIONS'
3810 ! include 'COMMON.LOCAL'
3811 ! include 'COMMON.VAR'
3812 ! include 'COMMON.CHAIN'
3813 ! include 'COMMON.INTERACT'
3814 ! include 'COMMON.IOUNITS'
3815 ! include 'COMMON.GEO'
3816 ! include 'COMMON.NAMES'
3817 ! include 'COMMON.CONTROL'
3818 ! include 'COMMON.DISTFIT'
3819 ! include 'COMMON.SETUP'
3820 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3822 logical :: lprn=.true.,fail
3823 real(kind=8),dimension(3) :: e1,e2,e3
3824 real(kind=8) :: dcj,efree_temp
3825 character(len=3) :: seq,res,res2
3826 character(len=5) :: atom
3827 character(len=80) :: card
3828 real(kind=8),dimension(3,40) :: sccor
3829 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3830 integer :: isugar,molecprev,firstion
3831 character*1 :: sugar
3833 real(kind=8),dimension(3) :: ccc
3835 integer,dimension(2,maxres/3) :: hfrag_alloc
3836 integer,dimension(4,maxres/3) :: bfrag_alloc
3837 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3838 real(kind=8),dimension(:,:), allocatable :: c_temporary
3839 integer,dimension(:,:) , allocatable :: itype_temporary
3840 integer,dimension(:),allocatable :: istype_temp
3847 ! write (2,*) "UNRES_PDB",unres_pdb
3867 !-----------------------------
3868 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3869 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3870 if(.not. allocated(istype)) allocate(istype(maxres))
3872 read (ipdbin,'(a80)',end=10) card
3873 write (iout,'(a)') card
3874 if (card(:5).eq.'HELIX') then
3877 read(card(22:25),*) hfrag(1,nhfrag)
3878 read(card(34:37),*) hfrag(2,nhfrag)
3880 if (card(:5).eq.'SHEET') then
3883 read(card(24:26),*) bfrag(1,nbfrag)
3884 read(card(35:37),*) bfrag(2,nbfrag)
3885 !rc----------------------------------------
3886 !rc to be corrected !!!
3887 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3888 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3889 !rc----------------------------------------
3891 if (card(:3).eq.'END') then
3893 else if (card(:3).eq.'TER') then
3898 itype(ires_old,molecule)=ntyp1_molec(molecule)
3899 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3900 nres_molec(molecule)=nres_molec(molecule)+2
3901 ! if (molecule.eq.4) ires=ires+2
3903 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3906 dc(j,ires)=sccor(j,iii)
3909 call sccenter(ires,iii,sccor)
3913 else if (card(:3).eq.'BRA') then
3917 itype(ires,molecule)=ntyp1_molec(molecule)-1
3918 nres_molec(molecule)=nres_molec(molecule)+1
3922 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3923 ! Fish out the ATOM cards.
3924 ! write(iout,*) 'card',card(1:20)
3925 ! print *,"ATU ",card(1:6), CARD(3:6)
3927 if (index(card(1:4),'ATOM').gt.0) then
3928 read (card(12:16),*) atom
3929 ! write (iout,*) "! ",atom," !",ires
3930 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3931 if (card(14:16).eq.'LIP') then
3936 nres_molec(molecule)=nres_molec(molecule)+1
3937 itype(ires,molecule)=ntyp1_molec(molecule)
3946 nres_molec(molecule)=nres_molec(molecule)+1
3947 read (card(18:20),'(a3)') res
3949 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3950 if (UNRES_PDB) then!
3951 itype(ires,molecule)=rescode(ires,res,0,molecule)
3953 itype(ires,molecule)=rescode_lip(res,ires)
3956 read (card(23:26),*) ires
3957 read (card(18:20),'(a3)') res
3958 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3959 ! & " ires_old",ires_old
3960 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3961 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3962 if (ires-ishift+ishift1.ne.ires_old) then
3963 ! Calculate the CM of the preceding residue.
3964 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3966 ! write (iout,*) "Calculating sidechain center iii",iii
3969 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3972 call sccenter(ires_old,iii,sccor)
3976 ! Start new residue.
3977 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3980 else if (ibeg.eq.1) then
3981 write (iout,*) "BEG ires",ires
3983 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3986 nres_molec(molecule)=nres_molec(molecule)+1
3988 ires=ires-ishift+ishift1
3990 ! write (iout,*) "ishift",ishift," ires",ires,&
3991 ! " ires_old",ires_old
3993 else if (ibeg.eq.2) then
3995 ishift=-ires_old+ires-1 !!!!!
3996 ishift1=ishift1-1 !!!!!
3997 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3998 ires=ires-ishift+ishift1
3999 ! print *,ires,ishift,ishift1
4003 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
4004 ires=ires-ishift+ishift1
4007 ! print *,'atom',ires,atom
4008 if (res.eq.'ACE' .or. res.eq.'NHE') then
4011 if (atom.eq.'CA '.or.atom.eq.'N ') then
4013 itype(ires,molecule)=rescode(ires,res,0,molecule)
4015 ! nres_molec(molecule)=nres_molec(molecule)+1
4019 itype(ires,molecule)=rescode(ires,res2,0,molecule)
4020 ! nres_molec(molecule)=nres_molec(molecule)+1
4021 read (card(19:19),'(a1)') sugar
4022 isugar=sugarcode(sugar,ires)
4023 ! if (ibeg.eq.1) then
4027 ! print *,"ires=",ires,istype(ires)
4033 ires=ires-ishift+ishift1
4035 ! write (iout,*) "ires_old",ires_old," ires",ires
4036 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
4039 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
4040 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
4041 res.eq.'NHE'.and.atom(:2).eq.'HN') then
4042 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4043 ! print *,ires,ishift,ishift1
4044 ! write (iout,*) "backbone ",atom
4046 write (iout,'(2i3,2x,a,3f8.3)') &
4047 ires,itype(ires,1),res,(c(j,ires),j=1,3)
4050 nres_molec(molecule)=nres_molec(molecule)+1
4052 sccor(j,iii)=c(j,ires)
4054 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
4055 atom.eq."C2'" .or. atom.eq."C3'" &
4056 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
4057 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
4058 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
4059 ! print *,ires,ishift,ishift1
4063 ! sccor(j,iii)=c(j,ires)
4066 c(j,ires)=c(j,ires)+ccc(j)/5.0
4068 print *,counter,molecule
4069 if (counter.eq.5) then
4071 nres_molec(molecule)=nres_molec(molecule)+1
4074 ! sccor(j,iii)=c(j,ires)
4078 ! print *, "ATOM",atom(1:3)
4079 else if (atom.eq."C5'") then
4080 read (card(19:19),'(a1)') sugar
4081 isugar=sugarcode(sugar,ires)
4086 ! print *,ires,istype(ires)
4090 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4091 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4092 nres_molec(molecule)=nres_molec(molecule)+1
4093 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4097 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4099 else if ((atom.eq."C1'").and.unres_pdb) then
4101 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4102 ! write (*,*) card(23:27),ires,itype(ires,1)
4103 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4104 atom.ne.'N' .and. atom.ne.'C' .and. &
4105 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4106 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4107 .and. atom.ne.'P '.and. &
4108 atom(1:1).ne.'H' .and. &
4109 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4111 ! write (iout,*) "sidechain ",atom
4112 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4113 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4114 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4116 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4120 ! print *,"IONS",ions,card(1:6)
4121 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4122 if (firstion.eq.0) then
4126 dc(j,ires)=sccor(j,iii)
4129 call sccenter(ires,iii,sccor)
4132 read (card(12:16),*) atom
4133 print *,"HETATOM", atom(1:2)
4134 read (card(18:20),'(a3)') res
4135 if (atom(3:3).eq.'H') cycle
4136 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4137 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4138 .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
4139 (atom(1:2).eq.'O ')) then
4141 print *,"I have water"
4142 if (molecule.ne.5) molecprev=molecule
4144 nres_molec(molecule)=nres_molec(molecule)+1
4145 print *,"HERE",nres_molec(molecule)
4146 if (res.ne.'WAT') res=res(2:3)//' '
4147 itype(ires,molecule)=rescode(ires,res,0,molecule)
4148 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4152 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4153 if (ires.eq.0) return
4154 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4157 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4159 nres_molec(molecule)=nres_molec(molecule)-2
4160 print *,'I have',nres, nres_molec(:)
4162 do k=1,4 ! ions are without dummy
4163 if (nres_molec(k).eq.0) cycle
4165 ! write (iout,*) i,itype(i,1)
4166 ! if (itype(i,1).eq.ntyp1) then
4167 ! write (iout,*) "dummy",i,itype(i,1)
4169 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4170 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4174 if (itype(i,k).eq.ntyp1_molec(k)) then
4175 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4176 if (itype(i+2,k).eq.0) then
4177 ! print *,"masz sieczke"
4179 if (itype(i+2,j).ne.0) then
4181 itype(i+1,j)=ntyp1_molec(j)
4182 nres_molec(k)=nres_molec(k)-1
4183 nres_molec(j)=nres_molec(j)+1
4189 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4190 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4191 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4192 ! if (unres_pdb) then
4193 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4194 ! print *,i,'tu dochodze'
4195 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4203 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4207 dcj=(c(j,i-2)-c(j,i-3))/2.0
4208 if (dcj.eq.0) dcj=1.23591524223
4213 else !itype(i+1,1).eq.ntyp1
4214 ! if (unres_pdb) then
4215 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4216 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4223 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4224 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4228 dcj=(c(j,i+3)-c(j,i+2))/2.0
4229 if (dcj.eq.0) dcj=1.23591524223
4234 endif !itype(i+1,1).eq.ntyp1
4235 endif !itype.eq.ntyp1
4239 ! Calculate the CM of the last side chain.
4243 dc(j,ires)=sccor(j,iii)
4246 call sccenter(ires,iii,sccor)
4252 ! print *,"molecule",molecule
4253 if ((itype(nres,1).ne.10)) then
4256 if (molecule.eq.5) molecule=molecprev
4257 itype(nres,molecule)=ntyp1_molec(molecule)
4258 nres_molec(molecule)=nres_molec(molecule)+1
4259 ! if (unres_pdb) then
4260 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4261 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4268 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4272 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4273 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4274 c(j,2*nres)=c(j,nres)
4278 ! print *,'I have',nres, nres_molec(:)
4280 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4282 if (nres.ne.nres0) then
4283 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4285 stop "Error nres value in WHAM input"
4288 !---------------------------------
4289 !el reallocate tables
4292 ! hfrag_alloc(j,i)=hfrag(j,i)
4295 ! bfrag_alloc(j,i)=bfrag(j,i)
4301 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4302 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4303 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4304 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4308 ! hfrag(j,i)=hfrag_alloc(j,i)
4313 ! bfrag(j,i)=bfrag_alloc(j,i)
4316 !el end reallocate tables
4317 !---------------------------------
4325 c(j,2*nres)=c(j,nres)
4328 if (itype(1,1).eq.ntyp1) then
4332 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4333 call refsys(2,3,4,e1,e2,e3,fail)
4340 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4341 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4345 dcj=(c(j,4)-c(j,3))/2.0
4351 ! First lets assign correct dummy to correct type of chain
4353 if (itype(1,1).eq.ntyp1) then
4354 if (itype(2,1).eq.0) then
4356 if (itype(2,j).ne.0) then
4358 itype(1,j)=ntyp1_molec(j)
4359 nres_molec(1)=nres_molec(1)-1
4360 nres_molec(j)=nres_molec(j)+1
4367 print *,'I have',nres, nres_molec(:)
4369 ! Copy the coordinates to reference coordinates
4375 ! Calculate internal coordinates.
4377 write (iout,'(/a)') &
4378 "Cartesian coordinates of the reference structure"
4379 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4380 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4382 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4383 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4384 (c(j,ires+nres),j=1,3)
4387 ! znamy już nres wiec mozna alokowac tablice
4388 ! Calculate internal coordinates.
4389 if(me.eq.king.or..not.out1file)then
4390 write (iout,'(a)') &
4391 "Backbone and SC coordinates as read from the PDB"
4393 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4394 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4395 (c(j,nres+ires),j=1,3)
4398 ! NOW LETS ROCK! SORTING
4399 allocate(c_temporary(3,2*nres))
4400 allocate(itype_temporary(nres,5))
4401 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4402 if (.not.allocated(istype)) write(iout,*) &
4403 "SOMETHING WRONG WITH ISTYTPE"
4404 allocate(istype_temp(nres))
4405 itype_temporary(:,:)=0
4409 if (itype(i,k).ne.0) then
4411 c_temporary(j,seqalingbegin)=c(j,i)
4412 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4414 itype_temporary(seqalingbegin,k)=itype(i,k)
4415 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4416 istype_temp(seqalingbegin)=istype(i)
4417 molnum(seqalingbegin)=k
4418 seqalingbegin=seqalingbegin+1
4424 c(j,i)=c_temporary(j,i)
4426 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4432 itype(i,k)=itype_temporary(i,k)
4433 istype(i)=istype_temp(i)
4436 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4437 ! I have only ions now dummy atoms in the system
4439 itype(1,5)=itype(2,5)
4445 itype(i,5)=itype(i+1,5)
4452 nres_molec(1)=nres_molec(1)-1
4454 ! if (itype(1,1).eq.ntyp1) then
4457 ! if (unres_pdb) then
4458 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4459 ! call refsys(2,3,4,e1,e2,e3,fail)
4466 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4470 ! dcj=(c(j,4)-c(j,3))/2.0
4472 ! c(j,nres+1)=c(j,1)
4478 write (iout,'(/a)') &
4479 "Cartesian coordinates of the reference structure after sorting"
4480 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4481 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4483 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4484 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4485 (c(j,ires+nres),j=1,3)
4489 print *,seqalingbegin,nres
4490 if(.not.allocated(vbld)) then
4491 allocate(vbld(2*nres))
4496 if(.not.allocated(vbld_inv)) then
4497 allocate(vbld_inv(2*nres))
4503 if(.not.allocated(theta)) then
4504 allocate(theta(nres+2))
4508 if(.not.allocated(phi)) allocate(phi(nres+2))
4509 if(.not.allocated(alph)) allocate(alph(nres+2))
4510 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4511 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4512 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4513 if(.not.allocated(costtab)) allocate(costtab(nres))
4514 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4515 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4516 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4517 if(.not.allocated(xxref)) allocate(xxref(nres))
4518 if(.not.allocated(yyref)) allocate(yyref(nres))
4519 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4520 if(.not.allocated(dc_norm)) then
4521 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4522 allocate(dc_norm(3,0:2*nres+2))
4525 write(iout,*) "before int_from_cart"
4526 call int_from_cart(.true.,.false.)
4527 call sc_loc_geom(.false.)
4528 write(iout,*) "after int_from_cart"
4532 thetaref(i)=theta(i)
4535 write(iout,*) "after thetaref"
4543 dc(j,i)=c(j,i+1)-c(j,i)
4544 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4549 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4550 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4552 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4556 ! Copy the coordinates to reference coordinates
4557 ! Splits to single chain if occurs
4561 ! cref(j,i,cou)=c(j,i)
4567 ! chomo(j,i,k)=c(j,i)
4570 ! write(iout,*) "after chomo"
4572 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4573 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4574 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4575 !-----------------------------
4579 write (iout,*) "symetr", symetr
4582 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4584 ! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4585 ! chain_length=lll-1
4587 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4592 cref(j,i,cou)=c(j,i)
4593 cref(j,i+nres,cou)=c(j,i+nres)
4595 chain_rep(j,lll,kkk)=c(j,i)
4596 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4601 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4603 ! do kkk=1,chain_length
4604 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4608 ! makes copy of chains
4609 write (iout,*) "symetr", symetr
4614 ! if (symetr.gt.1) then
4615 ! call permut(symetr)
4621 ! write(iout,*) (tabperm(i,kkk),kkk=1,4)
4626 ! icha=tabperm(i,kkk)
4627 ! write (iout,*) i,icha
4628 ! do lll=1,chain_length
4630 ! if (cou.le.nres) then
4632 ! kupa=mod(lll,chain_length)
4633 ! iprzes=(kkk-1)*chain_length+lll
4634 ! if (kupa.eq.0) kupa=chain_length
4635 ! write (iout,*) "kupa", kupa
4636 ! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4637 ! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4644 !-koniec robienia kopii
4647 write (iout,*) "nowa struktura", nperm
4649 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4651 cref(3,i,kkk),cref(1,nres+i,kkk),&
4652 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4654 100 format (//' alpha-carbon coordinates ',&
4655 ' centroid coordinates'/ &
4656 ' ', 6X,'X',11X,'Y',11X,'Z', &
4657 10X,'X',11X,'Y',11X,'Z')
4658 110 format (a,'(',i5,')',6f12.5)
4664 bfrag(i,j)=bfrag(i,j)-ishift
4670 hfrag(i,j)=hfrag(i,j)-ishift
4676 end subroutine readpdb
4677 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4678 !-----------------------------------------------------------------------------
4680 !-----------------------------------------------------------------------------
4681 subroutine read_control
4695 use random, only: random_init
4696 ! implicit real*8 (a-h,o-z)
4697 ! include 'DIMENSIONS'
4699 use prng, only:prng_restart
4701 logical :: OKRandom!, prng_restart
4704 ! include 'COMMON.IOUNITS'
4705 ! include 'COMMON.TIME1'
4706 ! include 'COMMON.THREAD'
4707 ! include 'COMMON.SBRIDGE'
4708 ! include 'COMMON.CONTROL'
4709 ! include 'COMMON.MCM'
4710 ! include 'COMMON.MAP'
4711 ! include 'COMMON.HEADER'
4712 ! include 'COMMON.CSA'
4713 ! include 'COMMON.CHAIN'
4714 ! include 'COMMON.MUCA'
4715 ! include 'COMMON.MD'
4716 ! include 'COMMON.FFIELD'
4717 ! include 'COMMON.INTERACT'
4718 ! include 'COMMON.SETUP'
4719 !el integer :: KDIAG,ICORFL,IXDR
4720 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4721 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4722 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4723 ! character(len=80) :: ucase
4724 character(len=640) :: controlcard
4726 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4728 usampl=.false. ! the default value of usample should be 0
4729 ! write(iout,*) "KURWA2",usampl
4733 read (INP,'(a)') titel
4734 call card_concat(controlcard,.true.)
4735 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4736 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4737 call reada(controlcard,'SEED',seed,0.0D0)
4738 call random_init(seed)
4739 ! Set up the time limit (caution! The time must be input in minutes!)
4740 read_cart=index(controlcard,'READ_CART').gt.0
4741 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4742 call readi(controlcard,'SYM',symetr,1)
4743 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4744 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4745 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4746 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4747 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4748 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4749 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4750 call reada(controlcard,'DRMS',drms,0.1D0)
4751 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4752 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4753 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4754 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4755 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4756 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4757 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4758 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4759 write (iout,'(a,f10.1)')'DRMS = ',drms
4760 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4761 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4763 call readi(controlcard,'NZ_START',nz_start,0)
4764 call readi(controlcard,'NZ_END',nz_end,0)
4765 call readi(controlcard,'IZ_SC',iz_sc,0)
4766 timlim=60.0D0*timlim
4767 safety = 60.0d0*safety
4770 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4771 !C SHIELD keyword sets if the shielding effect of side-chains is used
4772 !C 0 denots no shielding is used all peptide are equally despite the
4773 !C solvent accesible area
4774 !C 1 the newly introduced function
4775 !C 2 reseved for further possible developement
4776 call readi(controlcard,'SHIELD',shield_mode,0)
4777 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4778 write(iout,*) "shield_mode",shield_mode
4779 call readi(controlcard,'TORMODE',tor_mode,0)
4780 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4781 write(iout,*) "torsional and valence angle mode",tor_mode
4783 !C Varibles set size of box
4784 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4785 protein=index(controlcard,"PROTEIN").gt.0
4786 ions=index(controlcard,"IONS").gt.0
4787 fodson=index(controlcard,"FODSON").gt.0
4788 call readi(controlcard,'OLDION',oldion,1)
4789 nucleic=index(controlcard,"NUCLEIC").gt.0
4790 write (iout,*) "with_theta_constr ",with_theta_constr
4791 AFMlog=(index(controlcard,'AFM'))
4792 selfguide=(index(controlcard,'SELFGUIDE'))
4793 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4794 call readi(controlcard,'GENCONSTR',genconstr,0)
4795 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4796 call reada(controlcard,'BOXY',boxysize,100.0d0)
4797 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4798 call readi(controlcard,'TUBEMOD',tubemode,0)
4799 print *,"SCELE",scelemode
4800 call readi(controlcard,"SCELEMODE",scelemode,0)
4801 print *,"SCELE",scelemode
4803 ! elemode = 0 is orignal UNRES electrostatics
4804 ! elemode = 1 is "Momo" potentials in progress
4805 ! elemode = 2 is in development EVALD
4808 write (iout,*) TUBEmode,"TUBEMODE"
4809 if (TUBEmode.gt.0) then
4810 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4811 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4812 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4813 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4814 call reada(controlcard,"VNANO",velnanoconst,0.0d0)
4815 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4816 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4817 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4818 buftubebot=bordtubebot+tubebufthick
4819 buftubetop=bordtubetop-tubebufthick
4822 ! CUTOFFF ON ELECTROSTATICS
4823 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4824 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4825 write(iout,*) "R_CUT_ELE=",r_cut_ele
4826 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4827 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4828 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4830 ! Lipidec parameters
4831 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4832 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4833 if (lipthick.gt.0.0d0) then
4834 bordliptop=(boxzsize+lipthick)/2.0
4835 bordlipbot=bordliptop-lipthick
4836 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4837 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4838 buflipbot=bordlipbot+lipbufthick
4839 bufliptop=bordliptop-lipbufthick
4840 if ((lipbufthick*2.0d0).gt.lipthick) &
4841 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4842 endif !lipthick.gt.0
4843 write(iout,*) "bordliptop=",bordliptop
4844 write(iout,*) "bordlipbot=",bordlipbot
4845 write(iout,*) "bufliptop=",bufliptop
4846 write(iout,*) "buflipbot=",buflipbot
4847 write (iout,*) "SHIELD MODE",shield_mode
4849 !C-------------------------
4850 minim=(index(controlcard,'MINIMIZE').gt.0)
4851 dccart=(index(controlcard,'CART').gt.0)
4852 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4853 overlapsc=.not.overlapsc
4854 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4855 searchsc=.not.searchsc
4856 sideadd=(index(controlcard,'SIDEADD').gt.0)
4857 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4858 outpdb=(index(controlcard,'PDBOUT').gt.0)
4859 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4860 pdbref=(index(controlcard,'PDBREF').gt.0)
4861 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4862 indpdb=index(controlcard,'PDBSTART')
4863 extconf=(index(controlcard,'EXTCONF').gt.0)
4864 call readi(controlcard,'IPRINT',iprint,0)
4865 call readi(controlcard,'MAXGEN',maxgen,10000)
4866 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4867 call readi(controlcard,"KDIAG",kdiag,0)
4868 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4869 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4870 write (iout,*) "RESCALE_MODE",rescale_mode
4871 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4872 if (index(controlcard,'REGULAR').gt.0.0D0) then
4873 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4877 if (index(controlcard,'CHECKGRAD').gt.0) then
4879 if (index(controlcard,'CART').gt.0) then
4881 elseif (index(controlcard,'CARINT').gt.0) then
4886 elseif (index(controlcard,'THREAD').gt.0) then
4888 call readi(controlcard,'THREAD',nthread,0)
4889 if (nthread.gt.0) then
4890 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4893 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4894 stop 'Error termination in Read_Control.'
4896 else if (index(controlcard,'MCMA').gt.0) then
4898 else if (index(controlcard,'MCEE').gt.0) then
4900 else if (index(controlcard,'MULTCONF').gt.0) then
4902 else if (index(controlcard,'MAP').gt.0) then
4904 call readi(controlcard,'MAP',nmap,0)
4905 else if (index(controlcard,'CSA').gt.0) then
4907 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4909 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4912 !fcm else if (index(controlcard,'MCMF').gt.0) then
4914 else if (index(controlcard,'SOFTREG').gt.0) then
4916 else if (index(controlcard,'CHECK_BOND').gt.0) then
4918 else if (index(controlcard,'TEST').gt.0) then
4920 else if (index(controlcard,'MD').gt.0) then
4922 else if (index(controlcard,'RE ').gt.0) then
4926 lmuca=index(controlcard,'MUCA').gt.0
4927 call readi(controlcard,'MUCADYN',mucadyn,0)
4928 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4929 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4931 write (iout,*) 'MUCADYN=',mucadyn
4932 write (iout,*) 'MUCASMOOTH=',muca_smooth
4935 iscode=index(controlcard,'ONE_LETTER')
4936 indphi=index(controlcard,'PHI')
4937 indback=index(controlcard,'BACK')
4938 iranconf=index(controlcard,'RAND_CONF')
4939 i2ndstr=index(controlcard,'USE_SEC_PRED')
4940 gradout=index(controlcard,'GRADOUT').gt.0
4941 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4942 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4943 if (me.eq.king .or. .not.out1file ) &
4944 write (iout,*) "DISTCHAINMAX",distchainmax
4946 if(me.eq.king.or..not.out1file) &
4947 write (iout,'(2a)') diagmeth(kdiag),&
4948 ' routine used to diagonalize matrices.'
4949 if (shield_mode.gt.0) then
4950 pi=4.0D0*datan(1.0D0)
4951 !C VSolvSphere the volume of solving sphere
4953 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4954 !C there will be no distinction between proline peptide group and normal peptide
4955 !C group in case of shielding parameters
4956 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4957 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4958 write (iout,*) VSolvSphere,VSolvSphere_div
4959 !C long axis of side chain
4961 ! long_r_sidechain(i)=vbldsc0(1,i)
4962 ! short_r_sidechain(i)=sigma0(i)
4963 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4970 end subroutine read_control
4971 !-----------------------------------------------------------------------------
4972 subroutine read_REMDpar
4974 ! Read REMD settings
4981 use control_data, only:out1file
4982 ! implicit real*8 (a-h,o-z)
4983 ! include 'DIMENSIONS'
4984 ! include 'COMMON.IOUNITS'
4985 ! include 'COMMON.TIME1'
4986 ! include 'COMMON.MD'
4989 !el include 'COMMON.LANGEVIN'
4991 !el include 'COMMON.LANGEVIN.lang0'
4993 ! include 'COMMON.INTERACT'
4994 ! include 'COMMON.NAMES'
4995 ! include 'COMMON.GEO'
4996 ! include 'COMMON.REMD'
4997 ! include 'COMMON.CONTROL'
4998 ! include 'COMMON.SETUP'
4999 ! character(len=80) :: ucase
5000 character(len=320) :: controlcard
5001 character(len=3200) :: controlcard1
5002 integer :: iremd_m_total
5005 ! real(kind=8) :: var,ene
5007 if(me.eq.king.or..not.out1file) &
5008 write (iout,*) "REMD setup"
5010 call card_concat(controlcard,.true.)
5011 call readi(controlcard,"NREP",nrep,3)
5012 call readi(controlcard,"NSTEX",nstex,1000)
5013 call reada(controlcard,"RETMIN",retmin,10.0d0)
5014 call reada(controlcard,"RETMAX",retmax,1000.0d0)
5015 mremdsync=(index(controlcard,'SYNC').gt.0)
5016 call readi(controlcard,"NSYN",i_sync_step,100)
5017 restart1file=(index(controlcard,'REST1FILE').gt.0)
5018 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
5019 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
5020 if(max_cache_traj_use.gt.max_cache_traj) &
5021 max_cache_traj_use=max_cache_traj
5022 if(me.eq.king.or..not.out1file) then
5023 !d if (traj1file) then
5024 !rc caching is in testing - NTWX is not ignored
5025 !d write (iout,*) "NTWX value is ignored"
5026 !d write (iout,*) " trajectory is stored to one file by master"
5027 !d write (iout,*) " before exchange at NSTEX intervals"
5029 write (iout,*) "NREP= ",nrep
5030 write (iout,*) "NSTEX= ",nstex
5031 write (iout,*) "SYNC= ",mremdsync
5032 write (iout,*) "NSYN= ",i_sync_step
5033 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
5036 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
5037 if (index(controlcard,'TLIST').gt.0) then
5039 call card_concat(controlcard1,.true.)
5040 read(controlcard1,*) (remd_t(i),i=1,nrep)
5041 if(me.eq.king.or..not.out1file) &
5042 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
5045 if (index(controlcard,'MLIST').gt.0) then
5047 call card_concat(controlcard1,.true.)
5048 read(controlcard1,*) (remd_m(i),i=1,nrep)
5049 if(me.eq.king.or..not.out1file) then
5050 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
5053 iremd_m_total=iremd_m_total+remd_m(i)
5055 write (iout,*) 'Total number of replicas ',iremd_m_total
5058 if(me.eq.king.or..not.out1file) &
5059 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
5061 end subroutine read_REMDpar
5062 !-----------------------------------------------------------------------------
5063 subroutine read_MDpar
5067 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
5069 use geometry_data, only: pi
5071 ! implicit real*8 (a-h,o-z)
5072 ! include 'DIMENSIONS'
5073 ! include 'COMMON.IOUNITS'
5074 ! include 'COMMON.TIME1'
5075 ! include 'COMMON.MD'
5078 !el include 'COMMON.LANGEVIN'
5080 !el include 'COMMON.LANGEVIN.lang0'
5082 ! include 'COMMON.INTERACT'
5083 ! include 'COMMON.NAMES'
5084 ! include 'COMMON.GEO'
5085 ! include 'COMMON.SETUP'
5086 ! include 'COMMON.CONTROL'
5087 ! include 'COMMON.SPLITELE'
5088 ! character(len=80) :: ucase
5089 character(len=320) :: controlcard
5094 call card_concat(controlcard,.true.)
5095 call readi(controlcard,"NSTEP",n_timestep,1000000)
5096 call readi(controlcard,"NTWE",ntwe,100)
5097 call readi(controlcard,"NTWX",ntwx,1000)
5098 call readi(controlcard,"NFOD",nfodstep,100)
5099 call reada(controlcard,"DT",d_time,1.0d-1)
5100 call reada(controlcard,"DVMAX",dvmax,2.0d1)
5101 call reada(controlcard,"DAMAX",damax,1.0d1)
5102 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
5103 call readi(controlcard,"LANG",lang,0)
5104 RESPA = index(controlcard,"RESPA") .gt. 0
5105 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
5106 ntime_split0=ntime_split
5107 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5108 ntime_split0=ntime_split
5109 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5110 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5111 rest = index(controlcard,"REST").gt.0
5112 tbf = index(controlcard,"TBF").gt.0
5113 usampl = index(controlcard,"USAMPL").gt.0
5114 ! write(iout,*) "KURWA",usampl
5115 mdpdb = index(controlcard,"MDPDB").gt.0
5116 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5117 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5118 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5119 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5120 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5121 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5122 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5123 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5124 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5125 large = index(controlcard,"LARGE").gt.0
5126 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5127 rattle = index(controlcard,"RATTLE").gt.0
5128 preminim=(index(controlcard,'PREMINIM').gt.0)
5129 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5130 write (iout,*) "PREMINIM ",preminim
5131 dccart=(index(controlcard,'CART').gt.0)
5132 if (preminim) call read_minim
5133 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5139 if(me.eq.king.or..not.out1file) then
5141 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5143 write (iout,'(a)') "The units are:"
5144 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5145 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5146 " acceleration: angstrom/(48.9 fs)**2"
5147 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5149 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5150 write (iout,'(a60,f10.5,a)') &
5151 "Initial time step of numerical integration:",d_time,&
5153 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5155 write (iout,'(2a,i4,a)') &
5156 "A-MTS algorithm used; initial time step for fast-varying",&
5157 " short-range forces split into",ntime_split," steps."
5158 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5159 r_cut," lambda",rlamb
5161 write (iout,'(2a,f10.5)') &
5162 "Maximum acceleration threshold to reduce the time step",&
5163 "/increase split number:",damax
5164 write (iout,'(2a,f10.5)') &
5165 "Maximum predicted energy drift to reduce the timestep",&
5166 "/increase split number:",edriftmax
5167 write (iout,'(a60,f10.5)') &
5168 "Maximum velocity threshold to reduce velocities:",dvmax
5169 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5170 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5171 if (rattle) write (iout,'(a60)') &
5172 "Rattle algorithm used to constrain the virtual bonds"
5176 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5177 call reada(controlcard,"RWAT",rwat,1.4d0)
5178 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5179 surfarea=index(controlcard,"SURFAREA").gt.0
5180 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5181 if(me.eq.king.or..not.out1file)then
5182 write (iout,'(/a,$)') "Langevin dynamics calculation"
5184 write (iout,'(a/)') &
5185 " with direct integration of Langevin equations"
5186 else if (lang.eq.2) then
5187 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5188 else if (lang.eq.3) then
5189 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5190 else if (lang.eq.4) then
5191 write (iout,'(a/)') " in overdamped mode"
5193 write (iout,'(//a,i5)') &
5194 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5197 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5198 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5199 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5200 write (iout,'(a60,f10.5)') &
5201 "Scaling factor of the friction forces:",scal_fric
5202 if (surfarea) write (iout,'(2a,i10,a)') &
5203 "Friction coefficients will be scaled by solvent-accessible",&
5204 " surface area every",reset_fricmat," steps."
5206 ! Calculate friction coefficients and bounds of stochastic forces
5207 eta=6*pi*cPoise*etawat
5208 if(me.eq.king.or..not.out1file) &
5209 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5212 do j=1,5 !types of molecules
5213 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5214 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5216 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5217 do j=1,5 !types of molecules
5219 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5220 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5224 if(me.eq.king.or..not.out1file)then
5226 write (iout,'(/2a/)') &
5227 "Radii of site types and friction coefficients and std's of",&
5228 " stochastic forces of fully exposed sites"
5229 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5232 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5233 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5238 if(me.eq.king.or..not.out1file)then
5239 write (iout,'(a)') "Berendsen bath calculation"
5240 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5241 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5243 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5244 count_reset_moment," steps"
5246 write (iout,'(a,i10,a)') &
5247 "Velocities will be reset at random every",count_reset_vel,&
5251 if(me.eq.king.or..not.out1file) &
5252 write (iout,'(a31)') "Microcanonical mode calculation"
5254 if(me.eq.king.or..not.out1file)then
5255 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5257 write(iout,*) "MD running with constraints."
5258 write(iout,*) "Equilibration time ", eq_time, " mtus."
5259 write(iout,*) "Constraining ", nfrag," fragments."
5260 write(iout,*) "Length of each fragment, weight and q0:"
5262 write (iout,*) "Set of restraints #",iset
5264 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5265 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5267 write(iout,*) "constraints between ", npair, "fragments."
5268 write(iout,*) "constraint pairs, weights and q0:"
5270 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5271 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5273 write(iout,*) "angle constraints within ", nfrag_back,&
5274 "backbone fragments."
5275 write(iout,*) "fragment, weights:"
5277 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5278 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5279 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5282 iset=mod(kolor,nset)+1
5285 if(me.eq.king.or..not.out1file) &
5286 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5288 end subroutine read_MDpar
5289 !-----------------------------------------------------------------------------
5293 ! implicit real*8 (a-h,o-z)
5294 ! include 'DIMENSIONS'
5295 ! include 'COMMON.MAP'
5296 ! include 'COMMON.IOUNITS'
5297 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5298 character(len=80) :: mapcard !,ucase
5301 ! real(kind=8) :: var,ene
5304 read (inp,'(a)') mapcard
5305 mapcard=ucase(mapcard)
5306 if (index(mapcard,'PHI').gt.0) then
5308 else if (index(mapcard,'THE').gt.0) then
5310 else if (index(mapcard,'ALP').gt.0) then
5312 else if (index(mapcard,'OME').gt.0) then
5315 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5316 stop 'Error - illegal variable spec in MAP card.'
5318 call readi (mapcard,'RES1',res1(imap),0)
5319 call readi (mapcard,'RES2',res2(imap),0)
5320 if (res1(imap).eq.0) then
5321 res1(imap)=res2(imap)
5322 else if (res2(imap).eq.0) then
5323 res2(imap)=res1(imap)
5325 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5326 write (iout,'(a)') &
5327 'Error - illegal definition of variable group in MAP.'
5328 stop 'Error - illegal definition of variable group in MAP.'
5330 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5331 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5332 call readi(mapcard,'NSTEP',nstep(imap),0)
5333 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5334 write (iout,'(a)') &
5335 'Illegal boundary and/or step size specification in MAP.'
5336 stop 'Illegal boundary and/or step size specification in MAP.'
5340 end subroutine map_read
5341 !-----------------------------------------------------------------------------
5344 use control_data, only: vdisulf
5346 ! implicit real*8 (a-h,o-z)
5347 ! include 'DIMENSIONS'
5348 ! include 'COMMON.IOUNITS'
5349 ! include 'COMMON.GEO'
5350 ! include 'COMMON.CSA'
5351 ! include 'COMMON.BANK'
5352 ! include 'COMMON.CONTROL'
5353 ! character(len=80) :: ucase
5354 character(len=620) :: mcmcard
5356 ! integer :: ntf,ik,iw_pdb
5357 ! real(kind=8) :: var,ene
5359 call card_concat(mcmcard,.true.)
5361 call readi(mcmcard,'NCONF',nconf,50)
5362 call readi(mcmcard,'NADD',nadd,0)
5363 call readi(mcmcard,'JSTART',jstart,1)
5364 call readi(mcmcard,'JEND',jend,1)
5365 call readi(mcmcard,'NSTMAX',nstmax,500000)
5366 call readi(mcmcard,'N0',n0,1)
5367 call readi(mcmcard,'N1',n1,6)
5368 call readi(mcmcard,'N2',n2,4)
5369 call readi(mcmcard,'N3',n3,0)
5370 call readi(mcmcard,'N4',n4,0)
5371 call readi(mcmcard,'N5',n5,0)
5372 call readi(mcmcard,'N6',n6,10)
5373 call readi(mcmcard,'N7',n7,0)
5374 call readi(mcmcard,'N8',n8,0)
5375 call readi(mcmcard,'N9',n9,0)
5376 call readi(mcmcard,'N14',n14,0)
5377 call readi(mcmcard,'N15',n15,0)
5378 call readi(mcmcard,'N16',n16,0)
5379 call readi(mcmcard,'N17',n17,0)
5380 call readi(mcmcard,'N18',n18,0)
5382 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5384 call readi(mcmcard,'NDIFF',ndiff,2)
5385 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5386 call readi(mcmcard,'IS1',is1,1)
5387 call readi(mcmcard,'IS2',is2,8)
5388 call readi(mcmcard,'NRAN0',nran0,4)
5389 call readi(mcmcard,'NRAN1',nran1,2)
5390 call readi(mcmcard,'IRR',irr,1)
5391 call readi(mcmcard,'NSEED',nseed,20)
5392 call readi(mcmcard,'NTOTAL',ntotal,10000)
5393 call reada(mcmcard,'CUT1',cut1,2.0d0)
5394 call reada(mcmcard,'CUT2',cut2,5.0d0)
5395 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5396 call readi(mcmcard,'ICMAX',icmax,3)
5397 call readi(mcmcard,'IRESTART',irestart,0)
5398 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5401 call reada(mcmcard,'DELE',dele,20.0d0)
5402 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5403 call readi(mcmcard,'IREF',iref,0)
5404 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5405 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5406 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5407 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5408 write (iout,*) "NCONF_IN",nconf_in
5410 end subroutine csaread
5411 !-----------------------------------------------------------------------------
5415 use control_data, only: MaxMoveType
5418 ! implicit real*8 (a-h,o-z)
5419 ! include 'DIMENSIONS'
5420 ! include 'COMMON.MCM'
5421 ! include 'COMMON.MCE'
5422 ! include 'COMMON.IOUNITS'
5423 ! character(len=80) :: ucase
5424 character(len=320) :: mcmcard
5427 ! real(kind=8) :: var,ene
5429 call card_concat(mcmcard,.true.)
5430 call readi(mcmcard,'MAXACC',maxacc,100)
5431 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5432 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5433 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5434 call readi(mcmcard,'MAXREPM',maxrepm,200)
5435 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5436 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5437 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5438 call reada(mcmcard,'E_UP',e_up,5.0D0)
5439 call reada(mcmcard,'DELTE',delte,0.1D0)
5440 call readi(mcmcard,'NSWEEP',nsweep,5)
5441 call readi(mcmcard,'NSTEPH',nsteph,0)
5442 call readi(mcmcard,'NSTEPC',nstepc,0)
5443 call reada(mcmcard,'TMIN',tmin,298.0D0)
5444 call reada(mcmcard,'TMAX',tmax,298.0D0)
5445 call readi(mcmcard,'NWINDOW',nwindow,0)
5446 call readi(mcmcard,'PRINT_MC',print_mc,0)
5447 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5448 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5449 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5450 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5451 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5452 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5453 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5454 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5455 if (nwindow.gt.0) then
5456 allocate(winstart(nwindow)) !!el (maxres)
5457 allocate(winend(nwindow)) !!el
5458 allocate(winlen(nwindow)) !!el
5459 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5461 winlen(i)=winend(i)-winstart(i)+1
5464 if (tmax.lt.tmin) tmax=tmin
5465 if (tmax.eq.tmin) then
5469 if (nstepc.gt.0 .and. nsteph.gt.0) then
5470 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5471 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5473 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5474 ! Probabilities of different move types
5475 sumpro_type(0)=0.0D0
5476 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5477 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5478 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5479 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5480 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5481 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5482 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5484 print *,'i',i,' sumprotype',sumpro_type(i)
5485 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5486 print *,'i',i,' sumprotype',sumpro_type(i)
5489 end subroutine mcmread
5490 !-----------------------------------------------------------------------------
5491 subroutine read_minim
5494 ! implicit real*8 (a-h,o-z)
5495 ! include 'DIMENSIONS'
5496 ! include 'COMMON.MINIM'
5497 ! include 'COMMON.IOUNITS'
5498 ! character(len=80) :: ucase
5499 character(len=320) :: minimcard
5501 ! integer :: ntf,ik,iw_pdb
5502 ! real(kind=8) :: var,ene
5504 call card_concat(minimcard,.true.)
5505 call readi(minimcard,'MAXMIN',maxmin,2000)
5506 call readi(minimcard,'MAXFUN',maxfun,5000)
5507 call readi(minimcard,'MINMIN',minmin,maxmin)
5508 call readi(minimcard,'MINFUN',minfun,maxmin)
5509 call reada(minimcard,'TOLF',tolf,1.0D-2)
5510 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5511 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5512 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5513 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5514 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5515 'Options in energy minimization:'
5516 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5517 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5518 'MinMin:',MinMin,' MinFun:',MinFun,&
5519 ' TolF:',TolF,' RTolF:',RTolF
5521 end subroutine read_minim
5522 !-----------------------------------------------------------------------------
5523 subroutine openunits
5525 use MD_data, only: usampl
5528 use control_data, only:out1file
5529 use control, only: getenv_loc
5530 ! implicit real*8 (a-h,o-z)
5531 ! include 'DIMENSIONS'
5534 character(len=16) :: form,nodename
5535 integer :: nodelen,ierror,npos
5537 ! include 'COMMON.SETUP'
5538 ! include 'COMMON.IOUNITS'
5539 ! include 'COMMON.MD'
5540 ! include 'COMMON.CONTROL'
5541 integer :: lenpre,lenpot,lentmp !,ilen
5543 character(len=3) :: out1file_text !,ucase
5544 character(len=3) :: ll
5547 ! integer :: ntf,ik,iw_pdb
5548 ! real(kind=8) :: var,ene
5550 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5551 call getenv_loc("PREFIX",prefix)
5553 call getenv_loc("POT",pot)
5554 call getenv_loc("DIRTMP",tmpdir)
5555 call getenv_loc("CURDIR",curdir)
5556 call getenv_loc("OUT1FILE",out1file_text)
5557 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5558 out1file_text=ucase(out1file_text)
5559 if (out1file_text(1:1).eq."Y") then
5562 out1file=fg_rank.gt.0
5567 if (lentmp.gt.0) then
5568 write (*,'(80(1h!))')
5569 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5570 write (*,'(80(1h!))')
5571 write (*,*)"All output files will be on node /tmp directory."
5573 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5574 if (me.eq.king) then
5575 write (*,*) "The master node is ",nodename
5576 else if (fg_rank.eq.0) then
5577 write (*,*) "I am the CG slave node ",nodename
5579 write (*,*) "I am the FG slave node ",nodename
5582 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5583 lenpre = lentmp+lenpre+1
5585 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5586 ! Get the names and open the input files
5587 #if defined(WINIFL) || defined(WINPGI)
5588 open(1,file=pref_orig(:ilen(pref_orig))// &
5589 '.inp',status='old',readonly,shared)
5590 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5591 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5592 ! Get parameter filenames and open the parameter files.
5593 call getenv_loc('BONDPAR',bondname)
5594 open (ibond,file=bondname,status='old',readonly,shared)
5595 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5596 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5597 call getenv_loc('THETPAR',thetname)
5598 open (ithep,file=thetname,status='old',readonly,shared)
5599 call getenv_loc('ROTPAR',rotname)
5600 open (irotam,file=rotname,status='old',readonly,shared)
5601 call getenv_loc('TORPAR',torname)
5602 open (itorp,file=torname,status='old',readonly,shared)
5603 call getenv_loc('TORDPAR',tordname)
5604 open (itordp,file=tordname,status='old',readonly,shared)
5605 call getenv_loc('FOURIER',fouriername)
5606 open (ifourier,file=fouriername,status='old',readonly,shared)
5607 call getenv_loc('ELEPAR',elename)
5608 open (ielep,file=elename,status='old',readonly,shared)
5609 call getenv_loc('SIDEPAR',sidename)
5610 open (isidep,file=sidename,status='old',readonly,shared)
5612 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5613 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5614 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5615 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5616 call getenv_loc('TORPAR_NUCL',torname_nucl)
5617 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5618 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5619 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5620 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5621 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5624 #elif (defined CRAY) || (defined AIX)
5625 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5627 ! print *,"Processor",myrank," opened file 1"
5628 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5629 ! print *,"Processor",myrank," opened file 9"
5630 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5631 ! Get parameter filenames and open the parameter files.
5632 call getenv_loc('BONDPAR',bondname)
5633 open (ibond,file=bondname,status='old',action='read')
5634 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5635 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5637 ! print *,"Processor",myrank," opened file IBOND"
5638 call getenv_loc('THETPAR',thetname)
5639 open (ithep,file=thetname,status='old',action='read')
5640 ! print *,"Processor",myrank," opened file ITHEP"
5641 call getenv_loc('ROTPAR',rotname)
5642 open (irotam,file=rotname,status='old',action='read')
5644 call getenv_loc('ROTPAR_END',rotname_end)
5645 open (irotam_end,file=rotname_end,status='old',action='read')
5647 ! print *,"Processor",myrank," opened file IROTAM"
5648 call getenv_loc('TORPAR',torname)
5649 open (itorp,file=torname,status='old',action='read')
5650 ! print *,"Processor",myrank," opened file ITORP"
5651 call getenv_loc('TORDPAR',tordname)
5652 open (itordp,file=tordname,status='old',action='read')
5653 ! print *,"Processor",myrank," opened file ITORDP"
5654 call getenv_loc('SCCORPAR',sccorname)
5655 open (isccor,file=sccorname,status='old',action='read')
5656 ! print *,"Processor",myrank," opened file ISCCOR"
5657 call getenv_loc('FOURIER',fouriername)
5658 open (ifourier,file=fouriername,status='old',action='read')
5659 ! print *,"Processor",myrank," opened file IFOURIER"
5660 call getenv_loc('ELEPAR',elename)
5661 open (ielep,file=elename,status='old',action='read')
5662 ! print *,"Processor",myrank," opened file IELEP"
5663 call getenv_loc('SIDEPAR',sidename)
5664 open (isidep,file=sidename,status='old',action='read')
5666 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5667 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5668 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5669 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5670 call getenv_loc('TORPAR_NUCL',torname_nucl)
5671 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5672 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5673 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5674 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5675 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5677 call getenv_loc('LIPTRANPAR',liptranname)
5678 open (iliptranpar,file=liptranname,status='old',action='read')
5679 call getenv_loc('TUBEPAR',tubename)
5680 open (itube,file=tubename,status='old',action='read')
5681 call getenv_loc('IONPAR',ionname)
5682 open (iion,file=ionname,status='old',action='read')
5683 call getenv_loc('IONPAR_TRAN',iontranname)
5684 open (iiontran,file=iontranname,status='old',action='read')
5685 ! print *,"Processor",myrank," opened file ISIDEP"
5686 ! print *,"Processor",myrank," opened parameter files"
5688 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5689 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5690 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5691 ! Get parameter filenames and open the parameter files.
5692 call getenv_loc('BONDPAR',bondname)
5693 open (ibond,file=bondname,status='old')
5694 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5695 open (ibond_nucl,file=bondname_nucl,status='old')
5697 call getenv_loc('THETPAR',thetname)
5698 open (ithep,file=thetname,status='old')
5699 call getenv_loc('ROTPAR',rotname)
5700 open (irotam,file=rotname,status='old')
5702 call getenv_loc('ROTPAR_END',rotname_end)
5703 open (irotam_end,file=rotname_end,status='old')
5705 call getenv_loc('TORPAR',torname)
5706 open (itorp,file=torname,status='old')
5707 call getenv_loc('TORDPAR',tordname)
5708 open (itordp,file=tordname,status='old')
5709 call getenv_loc('SCCORPAR',sccorname)
5710 open (isccor,file=sccorname,status='old')
5711 call getenv_loc('FOURIER',fouriername)
5712 open (ifourier,file=fouriername,status='old')
5713 call getenv_loc('ELEPAR',elename)
5714 open (ielep,file=elename,status='old')
5715 call getenv_loc('SIDEPAR',sidename)
5716 open (isidep,file=sidename,status='old')
5718 open (ithep_nucl,file=thetname_nucl,status='old')
5719 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5720 open (irotam_nucl,file=rotname_nucl,status='old')
5721 call getenv_loc('TORPAR_NUCL',torname_nucl)
5722 open (itorp_nucl,file=torname_nucl,status='old')
5723 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5724 open (itordp_nucl,file=tordname_nucl,status='old')
5725 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5726 open (isidep_nucl,file=sidename_nucl,status='old')
5728 call getenv_loc('LIPTRANPAR',liptranname)
5729 open (iliptranpar,file=liptranname,status='old')
5730 call getenv_loc('TUBEPAR',tubename)
5731 open (itube,file=tubename,status='old')
5732 call getenv_loc('IONPAR',ionname)
5733 open (iion,file=ionname,status='old')
5734 call getenv_loc('IONPAR_NUCL',ionnuclname)
5735 open (iionnucl,file=ionnuclname,status='old')
5736 call getenv_loc('IONPAR_TRAN',iontranname)
5737 open (iiontran,file=iontranname,status='old')
5738 call getenv_loc('WATWAT',iwaterwatername)
5739 open (iwaterwater,file=iwaterwatername,status='old')
5740 call getenv_loc('WATPROT',iwaterscname)
5741 open (iwatersc,file=iwaterscname,status='old')
5744 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5746 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5747 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5748 ! Get parameter filenames and open the parameter files.
5749 call getenv_loc('BONDPAR',bondname)
5750 open (ibond,file=bondname,status='old',action='read')
5751 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5752 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5753 call getenv_loc('THETPAR',thetname)
5754 open (ithep,file=thetname,status='old',action='read')
5755 call getenv_loc('ROTPAR',rotname)
5756 open (irotam,file=rotname,status='old',action='read')
5758 call getenv_loc('ROTPAR_END',rotname_end)
5759 open (irotam_end,file=rotname_end,status='old',action='read')
5761 call getenv_loc('TORPAR',torname)
5762 open (itorp,file=torname,status='old',action='read')
5763 call getenv_loc('TORDPAR',tordname)
5764 open (itordp,file=tordname,status='old',action='read')
5765 call getenv_loc('SCCORPAR',sccorname)
5766 open (isccor,file=sccorname,status='old',action='read')
5768 call getenv_loc('THETPARPDB',thetname_pdb)
5769 print *,"thetname_pdb ",thetname_pdb
5770 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5771 print *,ithep_pdb," opened"
5773 call getenv_loc('FOURIER',fouriername)
5774 open (ifourier,file=fouriername,status='old',readonly)
5775 call getenv_loc('ELEPAR',elename)
5776 open (ielep,file=elename,status='old',readonly)
5777 call getenv_loc('SIDEPAR',sidename)
5778 open (isidep,file=sidename,status='old',readonly)
5780 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5781 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5782 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5783 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5784 call getenv_loc('TORPAR_NUCL',torname_nucl)
5785 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5786 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5787 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5788 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5789 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5790 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5791 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5792 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5793 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5794 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5795 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5796 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5797 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5800 call getenv_loc('LIPTRANPAR',liptranname)
5801 open (iliptranpar,file=liptranname,status='old',action='read')
5802 call getenv_loc('LIPBOND',lipbondname)
5803 open (ilipbond,file=lipbondname,status='old',action='read')
5804 call getenv_loc('LIPNONBOND',lipnonbondname)
5805 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5806 call getenv_loc('TUBEPAR',tubename)
5807 open (itube,file=tubename,status='old',action='read')
5808 call getenv_loc('IONPAR',ionname)
5809 open (iion,file=ionname,status='old',action='read')
5810 call getenv_loc('IONPAR_NUCL',ionnuclname)
5811 open (iionnucl,file=ionnuclname,status='old',action='read')
5812 call getenv_loc('IONPAR_TRAN',iontranname)
5813 open (iiontran,file=iontranname,status='old',action='read')
5814 call getenv_loc('WATWAT',iwaterwatername)
5815 open (iwaterwater,file=iwaterwatername,status='old',action='read')
5816 call getenv_loc('WATPROT',iwaterscname)
5817 open (iwatersc,file=iwaterscname,status='old',action='read')
5820 call getenv_loc('ROTPARPDB',rotname_pdb)
5821 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5824 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5825 #if defined(WINIFL) || defined(WINPGI)
5826 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5827 #elif (defined CRAY) || (defined AIX)
5828 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5830 open (iscpp_nucl,file=scpname_nucl,status='old')
5832 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5837 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5838 ! Use -DOLDSCP to use hard-coded constants instead.
5840 call getenv_loc('SCPPAR',scpname)
5841 #if defined(WINIFL) || defined(WINPGI)
5842 open (iscpp,file=scpname,status='old',readonly,shared)
5843 #elif (defined CRAY) || (defined AIX)
5844 open (iscpp,file=scpname,status='old',action='read')
5846 open (iscpp,file=scpname,status='old')
5848 open (iscpp,file=scpname,status='old',action='read')
5851 call getenv_loc('PATTERN',patname)
5852 #if defined(WINIFL) || defined(WINPGI)
5853 open (icbase,file=patname,status='old',readonly,shared)
5854 #elif (defined CRAY) || (defined AIX)
5855 open (icbase,file=patname,status='old',action='read')
5857 open (icbase,file=patname,status='old')
5859 open (icbase,file=patname,status='old',action='read')
5862 ! Open output file only for CG processes
5863 ! print *,"Processor",myrank," fg_rank",fg_rank
5864 if (fg_rank.eq.0) then
5866 if (nodes.eq.1) then
5869 npos = dlog10(dfloat(nodes-1))+1
5871 if (npos.lt.3) npos=3
5872 write (liczba,'(i1)') npos
5873 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5875 write (liczba,form) me
5876 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5877 liczba(:ilen(liczba))
5878 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5880 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5882 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5883 liczba(:ilen(liczba))//'.mol2'
5884 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5885 liczba(:ilen(liczba))//'.stat'
5887 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5888 //liczba(:ilen(liczba))//'.stat')
5889 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5892 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5893 liczba(:ilen(liczba))//'.const'
5898 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5899 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5900 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5901 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5902 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5904 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5906 rest2name=prefix(:ilen(prefix))//'.rst'
5908 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5911 #if defined(AIX) || defined(PGI)
5912 if (me.eq.king .or. .not. out1file) &
5913 open(iout,file=outname,status='unknown')
5915 if (fg_rank.gt.0) then
5916 write (liczba,'(i3.3)') myrank/nfgtasks
5917 write (ll,'(bz,i3.3)') fg_rank
5918 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5923 open(igeom,file=intname,status='unknown',position='append')
5924 open(ipdb,file=pdbname,status='unknown')
5925 open(imol2,file=mol2name,status='unknown')
5926 open(istat,file=statname,status='unknown',position='append')
5928 !1out open(iout,file=outname,status='unknown')
5931 if (me.eq.king .or. .not.out1file) &
5932 open(iout,file=outname,status='unknown')
5934 if (fg_rank.gt.0) then
5935 write (liczba,'(i3.3)') myrank/nfgtasks
5936 write (ll,'(bz,i3.3)') fg_rank
5937 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5942 open(igeom,file=intname,status='unknown',access='append')
5943 open(ipdb,file=pdbname,status='unknown')
5944 open(imol2,file=mol2name,status='unknown')
5945 open(istat,file=statname,status='unknown',access='append')
5947 !1out open(iout,file=outname,status='unknown')
5950 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5951 csa_seed=prefix(:lenpre)//'.CSA.seed'
5952 csa_history=prefix(:lenpre)//'.CSA.history'
5953 csa_bank=prefix(:lenpre)//'.CSA.bank'
5954 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5955 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5956 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5957 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5958 csa_int=prefix(:lenpre)//'.int'
5959 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5960 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5961 csa_in=prefix(:lenpre)//'.CSA.in'
5962 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5965 write (iout,'(80(1h-))')
5966 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5967 write (iout,'(80(1h-))')
5968 write (iout,*) "Input file : ",&
5969 pref_orig(:ilen(pref_orig))//'.inp'
5970 write (iout,*) "Output file : ",&
5971 outname(:ilen(outname))
5973 write (iout,*) "Sidechain potential file : ",&
5974 sidename(:ilen(sidename))
5976 write (iout,*) "SCp potential file : ",&
5977 scpname(:ilen(scpname))
5979 write (iout,*) "Electrostatic potential file : ",&
5980 elename(:ilen(elename))
5981 write (iout,*) "Cumulant coefficient file : ",&
5982 fouriername(:ilen(fouriername))
5983 write (iout,*) "Torsional parameter file : ",&
5984 torname(:ilen(torname))
5985 write (iout,*) "Double torsional parameter file : ",&
5986 tordname(:ilen(tordname))
5987 write (iout,*) "SCCOR parameter file : ",&
5988 sccorname(:ilen(sccorname))
5989 write (iout,*) "Bond & inertia constant file : ",&
5990 bondname(:ilen(bondname))
5991 write (iout,*) "Bending parameter file : ",&
5992 thetname(:ilen(thetname))
5993 write (iout,*) "Rotamer parameter file : ",&
5994 rotname(:ilen(rotname))
5997 write (iout,*) "Thetpdb parameter file : ",&
5998 thetname_pdb(:ilen(thetname_pdb))
6001 write (iout,*) "Threading database : ",&
6002 patname(:ilen(patname))
6004 write (iout,*)" DIRTMP : ",&
6006 write (iout,'(80(1h-))')
6009 end subroutine openunits
6010 !-----------------------------------------------------------------------------
6013 use geometry_data, only: nres,dc
6015 ! implicit real*8 (a-h,o-z)
6016 ! include 'DIMENSIONS'
6017 ! include 'COMMON.CHAIN'
6018 ! include 'COMMON.IOUNITS'
6019 ! include 'COMMON.MD'
6022 ! real(kind=8) :: var,ene
6024 open(irest2,file=rest2name,status='unknown')
6025 read(irest2,*) totT,EK,potE,totE,t_bath
6028 ! AL 4/17/17: Now reading d_t(0,:) too
6030 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
6033 ! AL 4/17/17: Now reading d_c(0,:) too
6035 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
6038 read (irest2,*) iset
6042 end subroutine readrst
6043 !-----------------------------------------------------------------------------
6044 subroutine read_fragments
6048 use control_data, only:out1file
6051 ! implicit real*8 (a-h,o-z)
6052 ! include 'DIMENSIONS'
6056 ! include 'COMMON.SETUP'
6057 ! include 'COMMON.CHAIN'
6058 ! include 'COMMON.IOUNITS'
6059 ! include 'COMMON.MD'
6060 ! include 'COMMON.CONTROL'
6063 ! real(kind=8) :: var,ene
6065 read(inp,*) nset,nfrag,npair,nfrag_back
6067 !el from module energy
6068 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
6069 if(.not.allocated(wfrag_back)) then
6070 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6071 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6073 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
6074 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
6076 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
6077 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
6080 if(me.eq.king.or..not.out1file) &
6081 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
6082 " nfrag_back",nfrag_back
6084 read(inp,*) mset(iset)
6086 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
6088 if(me.eq.king.or..not.out1file) &
6089 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
6090 ifrag(2,i,iset), qinfrag(i,iset)
6093 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
6095 if(me.eq.king.or..not.out1file) &
6096 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
6097 ipair(2,i,iset), qinpair(i,iset)
6100 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6101 wfrag_back(3,i,iset),&
6102 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6103 if(me.eq.king.or..not.out1file) &
6104 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6105 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6109 end subroutine read_fragments
6110 !-----------------------------------------------------------------------------
6112 !-----------------------------------------------------------------------------
6116 ! implicit real*8 (a-h,o-z)
6117 ! include 'DIMENSIONS'
6118 ! include 'COMMON.CSA'
6119 ! include 'COMMON.BANK'
6120 ! include 'COMMON.IOUNITS'
6122 ! integer :: ntf,ik,iw_pdb
6123 ! real(kind=8) :: var,ene
6125 open(icsa_in,file=csa_in,status="old",err=100)
6126 read(icsa_in,*) nconf
6127 read(icsa_in,*) jstart,jend
6128 read(icsa_in,*) nstmax
6129 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6130 read(icsa_in,*) nran0,nran1,irr
6131 read(icsa_in,*) nseed
6132 read(icsa_in,*) ntotal,cut1,cut2
6133 read(icsa_in,*) estop
6134 read(icsa_in,*) icmax,irestart
6135 read(icsa_in,*) ntbankm,dele,difcut
6136 read(icsa_in,*) iref,rmscut,pnccut
6137 read(icsa_in,*) ndiff
6144 end subroutine csa_read
6145 !-----------------------------------------------------------------------------
6146 subroutine initial_write
6149 ! implicit real*8 (a-h,o-z)
6150 ! include 'DIMENSIONS'
6151 ! include 'COMMON.CSA'
6152 ! include 'COMMON.BANK'
6153 ! include 'COMMON.IOUNITS'
6155 ! integer :: ntf,ik,iw_pdb
6156 ! real(kind=8) :: var,ene
6158 open(icsa_seed,file=csa_seed,status="unknown")
6159 write(icsa_seed,*) "seed"
6161 #if defined(AIX) || defined(PGI)
6162 open(icsa_history,file=csa_history,status="unknown",&
6165 open(icsa_history,file=csa_history,status="unknown",&
6168 write(icsa_history,*) nconf
6169 write(icsa_history,*) jstart,jend
6170 write(icsa_history,*) nstmax
6171 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6172 write(icsa_history,*) nran0,nran1,irr
6173 write(icsa_history,*) nseed
6174 write(icsa_history,*) ntotal,cut1,cut2
6175 write(icsa_history,*) estop
6176 write(icsa_history,*) icmax,irestart
6177 write(icsa_history,*) ntbankm,dele,difcut
6178 write(icsa_history,*) iref,rmscut,pnccut
6179 write(icsa_history,*) ndiff
6181 write(icsa_history,*)
6184 open(icsa_bank1,file=csa_bank1,status="unknown")
6185 write(icsa_bank1,*) 0
6189 end subroutine initial_write
6190 !-----------------------------------------------------------------------------
6191 subroutine restart_write
6194 ! implicit real*8 (a-h,o-z)
6195 ! include 'DIMENSIONS'
6196 ! include 'COMMON.IOUNITS'
6197 ! include 'COMMON.CSA'
6198 ! include 'COMMON.BANK'
6200 ! integer :: ntf,ik,iw_pdb
6201 ! real(kind=8) :: var,ene
6203 #if defined(AIX) || defined(PGI)
6204 open(icsa_history,file=csa_history,position="append")
6206 open(icsa_history,file=csa_history,access="append")
6208 write(icsa_history,*)
6209 write(icsa_history,*) "This is restart"
6210 write(icsa_history,*)
6211 write(icsa_history,*) nconf
6212 write(icsa_history,*) jstart,jend
6213 write(icsa_history,*) nstmax
6214 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6215 write(icsa_history,*) nran0,nran1,irr
6216 write(icsa_history,*) nseed
6217 write(icsa_history,*) ntotal,cut1,cut2
6218 write(icsa_history,*) estop
6219 write(icsa_history,*) icmax,irestart
6220 write(icsa_history,*) ntbankm,dele,difcut
6221 write(icsa_history,*) iref,rmscut,pnccut
6222 write(icsa_history,*) ndiff
6223 write(icsa_history,*)
6224 write(icsa_history,*) "irestart is: ", irestart
6226 write(icsa_history,*)
6230 end subroutine restart_write
6231 !-----------------------------------------------------------------------------
6233 !-----------------------------------------------------------------------------
6234 subroutine write_pdb(npdb,titelloc,ee)
6236 ! implicit real*8 (a-h,o-z)
6237 ! include 'DIMENSIONS'
6238 ! include 'COMMON.IOUNITS'
6239 character(len=50) :: titelloc1
6240 character*(*) :: titelloc
6241 character(len=3) :: zahl
6242 character(len=5) :: liczba5
6244 integer :: npdb !,ilen
6248 ! real(kind=8) :: var,ene
6252 if (npdb.lt.1000) then
6253 call numstr(npdb,zahl)
6254 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6256 if (npdb.lt.10000) then
6257 write(liczba5,'(i1,i4)') 0,npdb
6259 write(liczba5,'(i5)') npdb
6261 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6263 call pdbout(ee,titelloc1,ipdb)
6266 end subroutine write_pdb
6267 !-----------------------------------------------------------------------------
6269 !-----------------------------------------------------------------------------
6270 subroutine write_thread_summary
6271 ! Thread the sequence through a database of known structures
6272 use control_data, only: refstr
6274 use energy_data, only: n_ene_comp
6276 ! implicit real*8 (a-h,o-z)
6277 ! include 'DIMENSIONS'
6279 use MPI_data !include 'COMMON.INFO'
6282 ! include 'COMMON.CONTROL'
6283 ! include 'COMMON.CHAIN'
6284 ! include 'COMMON.DBASE'
6285 ! include 'COMMON.INTERACT'
6286 ! include 'COMMON.VAR'
6287 ! include 'COMMON.THREAD'
6288 ! include 'COMMON.FFIELD'
6289 ! include 'COMMON.SBRIDGE'
6290 ! include 'COMMON.HEADER'
6291 ! include 'COMMON.NAMES'
6292 ! include 'COMMON.IOUNITS'
6293 ! include 'COMMON.TIME1'
6295 integer,dimension(maxthread) :: ip
6296 real(kind=8),dimension(0:n_ene) :: energia
6298 integer :: i,j,ii,jj,ipj,ik,kk,ist
6299 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6301 write (iout,'(30x,a/)') &
6302 ' *********** Summary threading statistics ************'
6303 write (iout,'(a)') 'Initial energies:'
6304 write (iout,'(a4,2x,a12,14a14,3a8)') &
6305 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6306 'RMSnat','NatCONT','NNCONT','RMS'
6307 ! Energy sort patterns
6312 enet=ener(n_ene-1,ip(i))
6315 if (ener(n_ene-1,ip(j)).lt.enet) then
6317 enet=ener(n_ene-1,ip(j))
6329 ist=nres_base(2,ii)+ipatt(2,i)
6331 energia(i)=ener0(kk,i)
6333 etot=ener0(n_ene_comp+1,i)
6334 rmsnat=ener0(n_ene_comp+2,i)
6335 rms=ener0(n_ene_comp+3,i)
6336 frac=ener0(n_ene_comp+4,i)
6337 frac_nn=ener0(n_ene_comp+5,i)
6340 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6341 i,str_nam(ii),ist+1,&
6342 (energia(print_order(kk)),kk=1,nprint_ene),&
6343 etot,rmsnat,frac,frac_nn,rms
6345 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6346 i,str_nam(ii),ist+1,&
6347 (energia(print_order(kk)),kk=1,nprint_ene),etot
6350 write (iout,'(//a)') 'Final energies:'
6351 write (iout,'(a4,2x,a12,17a14,3a8)') &
6352 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6353 'RMSnat','NatCONT','NNCONT','RMS'
6357 ist=nres_base(2,ii)+ipatt(2,i)
6359 energia(kk)=ener(kk,ik)
6361 etot=ener(n_ene_comp+1,i)
6362 rmsnat=ener(n_ene_comp+2,i)
6363 rms=ener(n_ene_comp+3,i)
6364 frac=ener(n_ene_comp+4,i)
6365 frac_nn=ener(n_ene_comp+5,i)
6366 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6367 i,str_nam(ii),ist+1,&
6368 (energia(print_order(kk)),kk=1,nprint_ene),&
6369 etot,rmsnat,frac,frac_nn,rms
6371 write (iout,'(/a/)') 'IEXAM array:'
6372 write (iout,'(i5)') nexcl
6374 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6376 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6377 'Max. time for threading step ',max_time_for_thread,&
6378 'Average time for threading step: ',ave_time_for_thread
6380 end subroutine write_thread_summary
6381 !-----------------------------------------------------------------------------
6382 subroutine write_stat_thread(ithread,ipattern,ist)
6384 use energy_data, only: n_ene_comp
6386 ! implicit real*8 (a-h,o-z)
6387 ! include "DIMENSIONS"
6388 ! include "COMMON.CONTROL"
6389 ! include "COMMON.IOUNITS"
6390 ! include "COMMON.THREAD"
6391 ! include "COMMON.FFIELD"
6392 ! include "COMMON.DBASE"
6393 ! include "COMMON.NAMES"
6394 real(kind=8),dimension(0:n_ene) :: energia
6396 integer :: ithread,ipattern,ist,i
6397 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6399 #if defined(AIX) || defined(PGI)
6400 open(istat,file=statname,position='append')
6402 open(istat,file=statname,access='append')
6405 energia(i)=ener(i,ithread)
6407 etot=ener(n_ene_comp+1,ithread)
6408 rmsnat=ener(n_ene_comp+2,ithread)
6409 rms=ener(n_ene_comp+3,ithread)
6410 frac=ener(n_ene_comp+4,ithread)
6411 frac_nn=ener(n_ene_comp+5,ithread)
6412 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6413 ithread,str_nam(ipattern),ist+1,&
6414 (energia(print_order(i)),i=1,nprint_ene),&
6415 etot,rmsnat,frac,frac_nn,rms
6418 end subroutine write_stat_thread
6419 !-----------------------------------------------------------------------------
6421 subroutine readpdb_template(k)
6422 ! Read the PDB file for read_constr_homology with read2sigma
6423 ! and convert the peptide geometry into virtual-chain geometry.
6425 ! include 'DIMENSIONS'
6426 ! include 'COMMON.LOCAL'
6427 ! include 'COMMON.VAR'
6428 ! include 'COMMON.CHAIN'
6429 ! include 'COMMON.INTERACT'
6430 ! include 'COMMON.IOUNITS'
6431 ! include 'COMMON.GEO'
6432 ! include 'COMMON.NAMES'
6433 ! include 'COMMON.CONTROL'
6434 ! include 'COMMON.FRAG'
6435 ! include 'COMMON.SETUP'
6436 use compare_data, only:nhfrag,nbfrag
6437 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6439 logical lprn /.false./,fail
6440 real(kind=8), dimension (3):: e1,e2,e3
6441 real(kind=8) :: dcj,efree_temp
6445 real(kind=8), dimension (3,40) :: sccor
6447 integer, dimension (:), allocatable :: iterter
6448 if(.not.allocated(iterter))allocate(iterter(nres))
6455 write (2,*) "UNRES_PDB",unres_pdb
6463 read (ipdbin,'(a80)',end=10) card
6464 if (card(:3).eq.'END') then
6466 else if (card(:3).eq.'TER') then
6469 itype(ires_old-1,1)=ntyp1
6470 iterter(ires_old-1)=1
6471 #if defined(WHAM_RUN) || defined(CLUSTER)
6472 if (ires_old.lt.nres) then
6474 itype(ires_old,1)=ntyp1
6476 #if defined(WHAM_RUN) || defined(CLUSTER)
6480 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6483 dc(j,ires)=sccor(j,iii)
6486 call sccenter(ires,iii,sccor)
6489 ! Fish out the ATOM cards.
6490 if (index(card(1:4),'ATOM').gt.0) then
6491 read (card(12:16),*) atom
6492 ! write (iout,*) "! ",atom," !",ires
6493 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6494 read (card(23:26),*) ires
6495 read (card(18:20),'(a3)') res
6496 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6497 ! & " ires_old",ires_old
6498 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6499 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6500 if (ires-ishift+ishift1.ne.ires_old) then
6501 ! Calculate the CM of the preceding residue.
6505 dc(j,ires_old)=sccor(j,iii)
6508 call sccenter(ires_old,iii,sccor)
6512 ! Start new residue.
6513 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6516 else if (ibeg.eq.1) then
6517 ! write (iout,*) "BEG ires",ires
6519 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6523 ires=ires-ishift+ishift1
6525 ! write (iout,*) "ishift",ishift," ires",ires,
6526 ! & " ires_old",ires_old
6527 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6529 else if (ibeg.eq.2) then
6531 ishift=-ires_old+ires-1
6533 ! write (iout,*) "New chain started",ires,ishift
6536 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6537 ires=ires-ishift+ishift1
6540 if (res.eq.'ACE' .or. res.eq.'NHE') then
6543 itype(ires,1)=rescode(ires,res,0,1)
6546 ires=ires-ishift+ishift1
6548 ! write (iout,*) "ires_old",ires_old," ires",ires
6549 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6552 ! write (2,*) "ires",ires," res ",res," ity",ity
6553 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6554 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6555 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6556 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6558 write (iout,'(2i3,2x,a,3f8.3)') &
6559 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6563 sccor(j,iii)=c(j,ires)
6565 if (ishift.ne.0) then
6566 ires_ca=ires+ishift-ishift1
6570 ! write (*,*) card(23:27),ires,itype(ires)
6571 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6572 atom.ne.'N' .and. atom.ne.'C' .and.&
6573 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6574 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6575 ! write (iout,*) "sidechain ",atom
6577 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6581 10 if(me.eq.king.or..not.out1file) &
6582 write (iout,'(a,i5)') ' Nres: ',ires
6583 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6586 write(2,*) "tutaj",ires,nres
6588 ! write (iout,*) i,itype(i),itype(i+1)
6589 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6590 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6591 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6592 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6593 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6595 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6596 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6603 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6607 dcj=(c(j,i-2)-c(j,i-3))/2.0
6608 if (dcj.eq.0) dcj=1.23591524223
6613 else !itype(i+1).eq.ntyp1
6615 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6616 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6623 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6627 dcj=(c(j,i+3)-c(j,i+2))/2.0
6628 if (dcj.eq.0) dcj=1.23591524223
6633 endif !itype(i+1).eq.ntyp1
6634 endif !itype.eq.ntyp1
6636 ! Calculate the CM of the last side chain.
6639 dc(j,ires)=sccor(j,iii)
6642 call sccenter(ires,iii,sccor)
6646 if (itype(nres,1).ne.10) then
6650 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6651 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6658 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6662 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6663 if (dcj.eq.0) dcj=1.23591524223
6664 c(j,nres)=c(j,nres-1)+dcj
6665 c(j,2*nres)=c(j,nres)
6676 c(j,2*nres)=c(j,nres)
6678 if (itype(1,1).eq.ntyp1) then
6682 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6683 call refsys(2,3,4,e1,e2,e3,fail)
6690 c(j,1)=c(j,2)-1.9d0*e2(j)
6694 dcj=(c(j,4)-c(j,3))/2.0
6700 ! Copy the coordinates to reference coordinates
6706 ! Calculate internal coordinates.
6707 if (out_template_coord) then
6708 write (iout,'(/a)') &
6709 "Cartesian coordinates of the reference structure"
6710 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6711 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6713 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6714 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6715 (c(j,ires+nres),j=1,3)
6718 ! Calculate internal coordinates.
6719 call int_from_cart(.true.,out_template_coord)
6720 call sc_loc_geom(.false.)
6722 thetaref(i)=theta(i)
6727 dc(j,i)=c(j,i+1)-c(j,i)
6728 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6733 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6734 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6736 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6737 ! & vbld_inv(i+nres)
6742 cref(j,i+nres,1)=c(j,i+nres)
6752 end subroutine readpdb_template
6753 !-----------------------------------------------------------------------------
6755 !-----------------------------------------------------------------------------
6756 end module io_config